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Abstract 

We describe a new method of calculation of generic multi-loop master integrals based on the 
numerical solution of systems of difference equations in one variable. We show algorithms for the 
construction of the systems using integration-by-parts identities and methods of solutions by means of 
expansions in factorial series and Laplace's transformation. We also describe new algorithms for the 
identification of master integrals and the reduction of generic Feynman integrals to master integrals, 
and procedures for generating and solving systems of differential equations in masses and momenta 
for master integrals. We apply our method to the calculation of the master integrals of massive 
vacuum and self-energy diagrams up to three loops and of massive vertex and box diagrams up to 
two loops. Implementation in a computer program of our approach is described. Important features 
of the implementation are: the ability to deal with hundreds of master integrals and the ability to 
obtain very high precision results expanded at will in the number of dimensions. 
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1 Introduction 



Nowadays the technique probably more popularly used for calculating the contribution of Feynman 
diagrams is that based on the integration- by-parts in D dimensions [1, 2]. After some algebraic oper- 
ations, like contracting Lorentz indices and calculating fermion traces, the contribution of the diagram 
is expressed as a combination of several Feynman integrals with different powers of numerators and 
denominators. 

This expression composed of many integrals is reduced to a combination of a limited number of 'master 
integrals' using recurrence relations obtained by combining the identities obtained by integration-by-parts. 
Then, the master integrals are calculated numerically or analytically, with some other method. 

As this technique is applied to more and more complicated diagrams, some difficulties appear. 

First, a working general algorithm for identifying master integrals and for obtaining such recurrence 
relations is not known at present^. Up to now, for each diagram laborious handwork was needed for 
obtaining such recurrence relations^. 

Second, the number of master integrals grows rapidly with the number of loops and legs of the 
diagrams. Considering for example the reduction to master integrals of the contribution to the g-2 
of the electron in QED, it is known that the one-, two-, and three- loop contributions are reduced by 
integration- by-parts to respectively, 1, 3 and 17 master integrals (sec [4, 5] for the analytical calculation 
of the three-loop contribution). At the four- loop level, to which integration- by-parts has been still not 
applied, one expects several hundreds of master integrals. 

Third, the calculation of a single multi-loop master integral is a difficult problem for which a general 
method applicable to any diagram, with any values of masses and momenta, able to provide high-precision 
values and suitable for automatic calculation, is not known at present. Up to now a variety of different 
methods has been used, according to the topology and the values of masses and momenta of the particular 
diagram considered. 

Therefore, to face problems like four-loop g-2, it arises the need of a completely automated approach to 
the c;alculation of Feynman integrals, and of finding out suitable new methods, algorithms and techniques 
of calculations which allow one to solve or avoid the above difficulties. 

In this paper we describe the methods and techniques to be used in such automated approach. In 
particular we present: 

1. A new method for determining the master integrals and for reducing generic Feynman integrals to 
master integrals, applicable to arbitrary Feynman diagrams. 

2. A new method for calculating master integrals based on the numerical solution of the recurrence 
relations provided by the integration-by-parts method, seen as linear difference equations in one 
index, applicable to arbitrary Feynman diagrams. 

3. A new method of generation and solution of systems of differential equations in masses and momenta 
for master integrals, applicable to arbitrary Feynman diagrams. 

The part more innovative and important of this work is the method of calculation of master integrals 
based on difference equations. As this mathematical topic appears (surprisingly) to be practically absent 
from the literature, in order to improve the intelligibility of the paper we will give an extensive discussion 
on this argument, including methods of solutions, techniques of calculations and examples of applications. 
Moreover, study of boundary conditions of difference equations will lead us to a detailed discussion of 
asymptotic behaviour of Feynman integrals for large powers of one denominator, another topic which has 
received very scarce attention in the literature. 

Another important point is the automation of the approach. All methods, algorithms and techniques 
developed in this work have been implemented in a comprehensive program called SYS, written on purpose, 
which, among other things, contains a simplified algebraic manipulator. The aim of this program is to 
calculate the value of a generic Feynman integral in completely automatic way, by reducing it to master 
integrals and by calculating the master integrals. The program turned out to be able to deal with 
diagrams with hundreds of master integrals and to obtain very high precision values (for example 100-200 
digits) expanded at will in e = (4 — D)/2, with all coefficients in numerical form, divergent terms included. 
Limitations of the current implementation will be described. 

^ An algorithm which in principle may solve this problem has been recently proposed in [3], but up to now no practical 
application was shown. 
^See section 2.3. 
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A crucial point is the test of the approach. By using SYS we have analyzed vacuum, self-energy, 
vertex and box diagrams up to three or two loops and we have calculated the values of master integrals 
for some values of masses and momenta, comparing them with already known results, when possible. We 
will show the results obtained. In particular we have calculated all the previously unknown single-scale 
massive three-loop self-energy master integrals. These numerical calculations, involving the calculation 
of hundreds of master integrals at the same time, allowed us to prove the reliability of the approach in 
real cases and to accumulate a considerable experience of calculation. 

The plan of the paper is as follows: In section 2 we illustrate an algorithm for the construction 
and solution of systems of integration-by-parts identities which allows one to reduce a generic Feynman 
integral to a combination of master integrals; this algorithm is the basis for the algorithms developed in 
the following sections. In section 3 we show how to obtain systems of difference equations in one variable 
satisfied by the master integrals. In section 4 we show how to find solutions of difi'erence equations using 
expansions in factorial series. In section 5 we discuss the methods to find the values of the arbitrary 
constants appearing in the solutions of the difference equations. In section 6 we describe the techniques 
used to sum the factorial series expansions. In section 7 we describe in detail the application of our 
methods to the one-loop self-energy diagram. In section 8 we illustrate an alternative method of solution 
of difference equations based on the Laplace's transformation. In section 9 we apply our methods to 
various diagrams from one to three loops. In section 10 we show how to deal with integrals with some 
particular values of masses and momenta, by solving systems of differential equations in masses and 
momenta. In section 11 we show some technical information on the computer program SYS used in the 
calculations. In section 12 we give our conclusions. 



2 Systems of identities between Feynman integrals 

After the introduction of the notation used for integrals and identities, in sections 2.3-2.4 we describe the 

new method used for solving the systems of intcgration-by-parts identities and for obtaining reduction 
to master integrals. The description is much detailed because the algorithm here shown is an essential 
part of this work, being the basis of all similar algorithms of solution of systems of identities which will 
be described in sections 3.2, 8.2 and 10.2. In section 2.5 we show how the algorithm can be used for 
determining the master integrals. In sections 2.6-2.7 we show how to use the algorithm for reducing inte- 
grals, and we describe its more interesting feature: the absence of intermediate integrals with increasing 
exponents of the denominators in the intermediate steps of the reduction to master integrals. 



2.1 Generalities 

Let us introduce some notations used throughout the paper. We consider a generic Feynman diagram with 
Nk loops, A^e external and Nd internal lines. The loop momenta are ki, i = 1, . . . ,Nk; the independent 
external momenta are Pi, i = 1, . . . , Np, where Np = Nf. — 1 (if A^e > 0) for total momentum conservation. 
The denominators of the propagators are Di = qf+mf, where qi is the momentum fiowing in the internal 
line i and nii is the mass of the internal line; each qi is a linear combination of the momenta {pj} and 
{kj}. The usual prescription qf +mf = q^ + — iO is implied if necessary. A generic Feynman integral 
regularized in f-dimensional euclidean momentum space has the form 

J[d''kr][d''k2]...[d''kN,]V^S , (1) 

where [d^ki] = d^ki/iT^/'^ and V~^s is the generic integrand 

Vjs = jv, , 7i > 0, 6iji > , (2) 

1 li=l i 

7 — {71, . . . y'jNd} ^-nd S = {Siji}. The numerator of Eq.(2) is a product of powers of all the possible 
scalar products involving the loop momenta k; the total number A^^p of such scalar products is 

N,p = NpNk + Nk{Nk + l)/2 . (3) 
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2.2 Algebraic and integration-by-parts identities 

Let us consider the generic integrand (2). For each denominator Dj we write the identity 

where {p ■ k)j indicates one scalar product involving loop momenta which appears in the expression of 
Dj, and Cj is its coefficient in the expression. The scalar products {p ■ k)^ must be chosen all different. 
These algebraic identities are applied in sequence to V^s and to the terms subsequently generated, more 
times if it is necessary. As a result, the original integrand V^s is transformed into a sum of new terms, 
each one not containing {p ■ k)- and Dj simultaneously, with the general form 

VLic.p= ' Jn r.a, , n<Na, aj,/3j>0, (5) 

llj=l ij 

where the subscript nia0 shows the dependence on the number n of denominators, their particular 
combination i = {ii, . . . ,?«}, and the exponents a = {ai, . . . ,a.n} and /? = ,f3N,p-n}- The 

symbols {p ■ k irred.)^-, j = 1,. . . , Nsp — n indicate the Nsp — n 'irreducible' scalar products which cannot 
be simplified further on by Eq.(4) with one of the denominators Di^, . . . , Di^ of V.^^^^] if n = Ngp the 
numerator of Eq.(5) is unity. We stress that one different set of irreducible scalar products corresponds 
to each different combination of denominators. 

Integrating by parts in D dimensions [1, 2] one can write the identities 

0, j = l,...,Nu, l = l,...,Np, 

(6) 

0, j,l = l,...,Nk , 

where V^.^^^^ is defined in Eq.(5). For each different V^.^^p Eq.(6) gives Nk{Np + Nk) identities. The 
ratio V^jQ,^ contains only irreducible scalar products (relative to the particular combination of denomi- 
nators); the calculation of the derivative and the contraction of the index /x form terms also containing 
reducible scalar products, which must be transformed into irreducible scalar products using the algebraic 
identities (4). As a final result, the identity will contain a linear combination of integrals of two kinds: 
integrals containing all the n denominators {Di^, Di^, . . . , Di^} appearing in and integrals with 

one denominator missing as effect of algebraic identities. 

2.3 A new method for solving the system of identities 

Each identity obtained from Eq.(6) is a linear combination of integrals like 

j[d''k,]...[d''kN,]VU^^, (7) 

with polynomials of degree zero or one in the number of dimensions D as coefficients, and where the 
integrand has the general form of Eq.(5). A generic set of these identities forms a homogeneous linear 
system of equations, with the integrals as unknowns. As it is well-known, the system is under-determined, 
and some integrals exist, the "master integrals" , whose values cannot be determined from the system. 

In the literature"^ it is usual to look for the general solution of the infinite system (6) in the form 
of combination of identities which, lowering and raising exponents, transform integrals of the form (1) 
or (7) into linear combinations of carefully chosen master integrals. Up to now, the methods used to 
find these identities have been based on laborious human analysis; as the number of these identities (and 
the difficulties encountered) grow rapidly with the number of denominators and loops, it becomes very 
difficult to analyze in this way diagrams beyond a certain limit. 

On the contrary, our approach is different, as it consists in the solution of systems made up of a finite 
number of identities. The identities (6) are generated explicitly using suitable V^'j^^^, with parameters n, 

^ For a review see [6]. In the program BUBBLES [7] reduction of systems of identities is performed by semi-automatic 
means. A different approach to solving integration-by-parts identities was developed in [8]. 



(4) 



3 



i, a, (3 taken from a large (but finite) set of values carefully chosen. The set of the generated identities 
forms a linear system with integrals as unknowns, which is solved using the well-known Gauss elimination 
method; the solution gives the expressions of the integrals as linear combinations of the master integrals 
with coefficients rational in D. The advantage of this approach is that it is simple, applicable without 
modifications to diagrams with any topology, and suitable for completely automatic calculations.^ It 
does not consist only in a mere "mechanical" approach because, as we will see in section 2.6, the solution 
of the whole system allows us to discover very useful identities of a kind a priori not expected; these 
identities allow one to avoid the appearing of intermediate integrals with increasing exponents of the 
denominators in the intermediate steps of the reduction to master integrals. 

Let lis now consider more in detail the solution of the system. The identities are generated and 
inserted in the system one at a time. Let X^^- CjWj = be an identity obtained from Eq.(6), where Wj 
are integrals and Cj the coefficients. The identities already existing in the system, expressing some of the 
integrals Wj in terms of other integrals are substituted in the new identity, which becomes '^'j^'j — 0- 
One particular integral is chosen between the integrals {W'j}, and the identity itself is rewritten as 
W/ = Ylj^i CjWj in order to express the particular integral in terms of the remaining integrals. Then the 
new identity is added to the system and the chosen integral W[ is substituted in the rest of the system. 

The choice of the integral W/ is carried out following an ordering^ of all the integrands V^^^^^ as 
functions of the parameters: the number n of denominators, the combination i of denominators, and 
the exponents. The ordering defines the priority of extraction between integrals. Priorities are arranged 
so that integrals with a higher number of denominators are extracted first, and expressed as linear 
combinations of the integrals with a lower number of denominators. Remaining details of ordering used 
are shown in the algorithm of the next section. The form of the master integrals depends on the choice 
of the ordering: see section 2.5. 

Now we consider the choice of the order of the generation and processing of the identities. The final 
solution of the system is independent of the order of processing, but the computing time is not. Each 
addition of a new identity to the system implies a substitution of an integral in all identities of the 
system which contains it; therefore the order must be carefully chosen in order to minimize the number 
of substitutions required. A bad choice may cause the computing time to blow up. 

A good choice of the ordering of the ratios V^^^^ appearing in Eq.(6) is the inverse of the above 
considered ordering of integrands used for the extraction of integrals. The identities corresponding to 
ratios V^^q,^ with the lowest priority of extraction (as integrands), therefore with the lowest number of 
denominators, are generated and processed first; then the identities with a higher number of denominators 
are processed. The last processed identities will have the highest number of denominators. This choice 
is efficient as long as the number n of denominators of the ratio V^j^^ used in Eq.(6) is equal to the 
number n' of denominators of the integral W[ extracted from the identity generated. In some cases 
(relatively rare, but of great importance, see section 2.6) one finds n' < n; the identities obtained with 
ratios V^n,^,^,, containing integrals with n' denominators, have already been inserted into the system, 
so that the "late" integral W[ must be substituted in a large number of identities, and that may require 
a considerable amount of time. 

It is useful to split the system of identities (6) into several subsystems. Each subsystem is made up 
of all the identities obtained by inserting in Eq.(6) terms V containing one particular combination i of 
denominators {Di-^, Di^, . . . , Di^}. The integrals with a number of denominators n less than the number 
of loops Nk of the diagram are null in the framework of the dimensional regularization, and must be not 
considered, so that the total number S of subsystems is 



A similar method, in very preliminary form, was developed by the author in [4] to reduce to master integrals the 
triple-cross vertex diagrams contribution to the g — 2 of the electron. In that work a system of about 100000 identities 
was built and solved, and this required some months of computer time; with our new method the same calculation may be 
performed in a fraction of hour. Some similar technique has been recently used in [9], to solve small systems of identities 
in order to obtain differential equations for two-loop massless box diagrams. 

^ The importance of a total ordering for the solution of systems of differential equations derived from integration-by-parts 
identities was remarked in [3]. 
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Let us now define the non- negative quantities Mp, Md, 



3=1 j=l 



(9) 



which are the total sum of the powers of the scalar products in the numerator and the sum of the 
powers of the denominators minus the number of denominators of the generic integrand V^j^^ of Eq.(5), 
respectively. Using this definition we can split the infinite set of the integrands VJ^j^^ into finite sets of 

. The set of the integrands 

. In the following sections 



integrands with equal n, Mp and that we indicate with the symbol 

with < Mp < a and < M^ < b will be indicated with the symbol [n; g '^ 
we will say in short that an integral (7) or an identity (6) belongs to a specified set, understanding that 
it is the corresponding V^j„^ which belongs to that set. 

2.4 The algorithm of solution of the system 

The algorithm of solution is the following: 
Algorithm 1 1. Let n = Nj.- 

2. Let ii = 1, i2 = 2, ... . in — n. 

3. Consider the combinations of n different denominators, chosen in the set {Di, Dn^}; let 
Di^ , Di^ , . . . , A„ be one of these combinations. 

4- Choose two integer non-negative constants Oi and bi, i = {ii, . . . ,i„}. 

5. Let Md = 0. 

6. Let Mp = 0. 

7. Consider the ratios of the kind of Eq.(5), containing the n denominators 
Di^,Di^, . . . , Di^ raised to some power; let W be one of these ratios: 



W{n, i, a, (3) 



11^=1 " • irred.)^^' 



where the non-negative exponents aj, j3j are such that 

Pj = Mp, E ("^- - 1) = , 



i.e., they are such that W belongs to the set 



n; 



Mp 



8. Set V^iQ,^ =W in Eq.(6) and generate all the integration-by-parts identities. 

9. Let j[d^ki] . . .[d,^kiq^] '^jCjWj = Q be one of the identities of Eq.(6), where Wj is a generic 
integrand with n or n — 1 denominators, then: 

(a) Substitute all the integrals already known in the left-hand side of the new identity; 

let /[d^fci] . . . [d^kNk] J2j c'jW'j = 6e the result. If the new identity is a linear combination 
of other identities of the system, go to step 10. 

(b) If the new identity is linearly independent, choose an integrand W/ from the identity. 



W/(n',i',a',/3') = ^^ 



Nsp-n' 



(p ■ k irred.)^,' 



belonging to 



M' 



, Nk < n' < n, following the order of priority: 
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i. the greatest number of the denominators n' ; 

a. the greatest M'^; 

in. the greatest M^; 

iv. the greatest i\, the greatest i'^, ... , the greatest i'^; 

V. the greatest a'l, the greatest a'^, ■ ■ ■ , the greatest a'^; 

vi. the greatest f3[, the greatest f3'2, ■ ■ ■ , the greatest P^^^-n'- 
(c) Substitute and add the following identity to the system: 

I [d^fci] . . . [d^fc^vj wi = - f [rf^fci] . . . [rf^fc^v.] • 
J J j^i 

10. Generate a new integration-by-parts identity among the Nh{Np + Nk) possible identities of Eq.(6) 
and go to step 9, otherwise continue. 



among 



11. Choose a new integrand W with different exponents a and (3, belonging to the set 
the (^"""'mp^''"^) ("^Md~^) e/emente of the set and go to step 8, otherwise continue. 

12. Mp = Mp + 1/ if Mp < Ui go to step 7. 

13. Md = Md + I; if Md < h go to step 6. 

14. Choose a new combination of indices ii < i2 < ■■ ■ < in among the (J^) possible combinations of 
numbers {1, 2, . . . , N^} and go to step 3, otherwise continue. 

15. n = n + 1; if n < N4 go to step 2, otherwise end. 

In this algorithm two arbitrary integer constants, ttj and 6,, must be chosen for each different combination 

of denominators {Di-^ , ■ ■ ■ , Di^}. These constants are cutoffs and define which identities must be included 
in the system and which must be excluded; limits the sum of the powers of scalar products in the 
numerator and bi limits the sum of the exponents of the denominators. 

With a suitable choice of the parameters and hi (see section 2.6), this algorithm allows one to 
reduce any given integral / to a sum of master integrals Bi 

L 

I = Y,riBi, (10) 
1=1 

where r; are rational functions of D, masses and scalar products of external momenta. 
2.5 Determination of master integrals 

Our algorithm of solution of the system of identities provides a general method for determining the master 
integrals: it suffices, for each different combination of n denominators {Di^, Di^, . . . , Di^} to build the 
subsystem with suitable and hi , to solve it using the algorithm, and to find the integrals belonging to 

n; q "^' which are not reduced to terms with n — 1 denominators. These should be master integrals. 
Of course, since we are dealing with a system made up with a finite number of identities, there is the 
possibility that a master integral is erroneously identified or not found because some essential identities 
have not been included in the system as they exceed the limits of generation. Explicit solutions of 
subsystems with increasing values of Oi and hi show that empirical 'thresholds' aj and hi exist, depending 
on the combination of denominators, such that all the subsystems built with > and hi > hi give a 
stable and correct identification of the master integrals; for the diagrams so far examined (see section 9) 
one finds < Si < 3 (it is equal to the maximum number of scalar products appearing in the master 
integrals) and bi = 0. Corresponding number of identities varies from some tens to some hundreds. 

The particular order of selection of the integrals used in step 9b of the algorithm favours for each 
different combination of denominators the choice of master integrals with all the denominators Dj raised 
to the first power and with the numerator containing increasing numbers of scalar products: first the 
scalar integral, with 1 as numerator, then integrals with one scalar product, and so on. 
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As result of the procedure just described we find the set of L master integrals Bi with the form 

where the combination of indices i and exponents j3 depends on the index I. The ordering Bi, B2, ... , 
Bl follows the ordering of integrals provided by the algorithm. 

2.6 Choice of constants and bii the "golden rule" 

The minimal values of and bi necessary to reduce one particular integral strongly depend on the 
structure of the integral and arc not easy to determine. Therefore wc limit ourselves to give some general 
rules, based on the experience acquired by using the algorithm, rather than rigorous results. Actually, 
these rules turned out to be very effective. 

Let us consider the system of identities obtained by inserting in Eq.(6) all the terms belonging 
to the set Gab 



Gab= U 
n=Nk 



0...a 



AT ... a 
. . . 



(12) 



containing all the possible S subsystems of identities (see Eq.(8)), and choosing cii = a and hi = b for all 
the combinations of denominators. The number of elements V^-^^^ forming the set Gab is given by 



( ^<^\ I — i^ + o\(n + b 



while the number of identities of the system is 

N^de (Gab) = Nk {Np + Nk) N,i,{Gab) . (14) 

Explicit calculations suggest the following empirical "golden rule" : the solution of the system provides 
identities which reduce any integral belonging to the set Gab to master integrals 

Integrals Gab > Master Integrals . (15) 

Identities Gab 

In order to better understand this rule, we have analyzed in some test cases the reduction of single 
integrals belonging to \_Nd] 'fc]' modifying the steps 1 and 15 of the algorithm so that the various 
subsystems of identities [n; q " "] are generated starting from Na denominators and decreasing n up to 
Nk denominators (note that the solution of the whole system or parts of it with such 'inverted' order is 
extremely time consuming) . For each value of n, intermediate expressions contain master integrals with a 

number of denominators from n to N^., and integrals belonging to the set n — 1; g° j,'^i , that is, integrals 

with one exponent of denominators increased by one compared with the original integral whichever is n. 
This is very different from what one obtains by performing the reduction one denominator at a time with 
recurrence relations which lower one index and raise another one, where the exponents of denominators 
of the integrals are observed to increase whenever the number of denominators decreases, causing the 
number of intermediate integrals to blow up. Our algorithm avoids this blowing-up. 

The quite unexpected behaviour of our algorithm seems to be due to the fact that solving the whole 
system we consider all the identities generated, not only those which reduce single integrals, but also 
those many "additional" identities apparently useless which reduce combinations of integrals with the 
same number of denominators. Their effect is to reduce systematically all the generated combinations 
of integrals with increased exponents of the denominators. The reason for this behaviour is not known; 
however, its systematic presence suggests that (hopefully) it may have some simple explanation. The 
presence of these "additional" identities is important, particularly if is large, because without them the 
number of identities to consider (and the number of intermediate integrals) would be orders of magnitude 
larger. 

The "golden rule" (15) turns out to be valid if a > ao and b > bo, where ao and 60 are some empirical 
'thresholds' depending on the structure of denominators of the diagram; for the diagrams examined in 
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section 9 wc found 1 < ao < 3 (values almost always equal to the maximum number of scalar products 
appearing in the master integrals) and bo = 0, 1. Values of the thresholds are unfortunately not known a 
priori; if the chosen values of a and b are below the threshold, the result of the solution of the system is 
that some integrals belonging to Gab arc not completely reduced to master integrals, and some non-master 
integrals with a few denominators still survive in the final expressions. In this case one may increase a or 
b, or alternatively one may suitably enlarge only the subsystems of identities corresponding to residual 
integrals. 

We found some very rare exceptions to the "golden rule" (15), when some denominators have zero 
mass; in these cases, whatever the values of a and b are, there are always few integrals belonging to 
[Nd] ^] (note, changing with the value of b) which are not completely reduced to master integrals. In the 
example of the next section we will show just one of these exceptions. 

Now, let us consider an important application: the reduction to master integrals of one combination 
of many integrals, for example, expressing the contribution of a diagram to some physical quantity. Let 
us suppose that all the integrals of this combination belong to a set Gg;, with some (minimal) a and b; 
if the values of a and b are over thresholds, according Eq.(15) the solution of the system will provide 
identities which reduce to master integrals all the integrals of the set. 



2.7 An example 

Let us consider the self-energy diagram with 5 denominators Di = (p — fci)^ + 1, D2 = {p — ki — k2)^ + 1, 
D3 = {p — ^2)^ + 1, D4 — kf, D5 = fe| and = —1. Following the notation of section 2.1, this diagram 
has Np = 1, Nk = 2, N4 = 5 and Nsp = 5. For instance, we want to transform the integral 

*1 (16) 



DiD2D3D4D5 



where [dk] = [d^ki] [d^k2], into a combination of master integrals. The pairs between scalar products 
and denominators used in the algebraic identities are (p- fci, Di), (/ci • A;2, D2), ip' k2, D3), {ki-ki,D4) and 
(^2 • k2, D5). We must identify the master integrals; therefore, following section 2.5, we build the system 
of identities with = 1 and bj = for all the combinations of denominators, and we look for the integrals 
which are not reduced. The master integrals turn out be: B123 = J[dk]/ D1D2D3, = J[dk]/ D3D4D5, 
B12 = J[dk]/D^D2, Bu = J[dk]/DiD3 and B23 = /[dfc]/£>2-D3. 

We consider the whole system made up of the identities of the set Gn, choosing a = b = 1; according 
to Eq.(14) it contains 1776 identities. Solving it with the algorithm 1 (by using the program SYS de- 
scribed in section 11) we find 1122 independent identities. Examining them, we find the identities which 
reduce to master integrals all the 291 integrals [2... 5; q"'\], integral J included, with the exception 
of the two integrals / [dk]/ D1D2D3DIDC, and / [dk]/ D1D2D3D4DI, whose reduction still contains the 
integrals / [dk]/D2DlDl^ and / [dk]/D3DlD-^; in order to complete the reduction of these two integrals, 
we may add the two sets of 30 identities [3; q "2] with these combinations of three denominators to the 
system. The remaining 831 identities reduce to master integrals complicated combinations of integrals 

[2...5;V]- . 

Wc want to illustrate how the mechanism of "additional" identities works. Let us consider the partial 
reduction of the integral J from 5 to 4 denominators. We consider the subsystem of the identities [5; q ° ^ ] , 
which is formed by 36 identities; solving it, we find: 

1. 6 identities reducing the integrals [5; to 4 denominators; 

2. 9 identities reducing combinations of integrals [5; 2] to 4 denominators; 

3. 20 identities between integrals [4; °] and [4; 2] with difi^erent combinations of denominators; 

The identities of the groups 2 and 3 are the "additional" ones. Among the identities of the group 1 we 
find 

J = J [dk] Qw^i + ^^^2) , (17) 
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We see that the original integral, which contains a square denominator (set [5; J]), is transformed into 
a combination of integrals with 4 denominators, with one or two square denominators (sets [4; Ij*] and 

Considering the reduction from 4 denominators to master integrals, if wc pciTisc the sohition of the 
subsystems [2 . . .4; q'"J] wc find 8 relations which transform each one of the terms of Wi into master 
integrals; on the contrary we do not find relations which transform single terms of W2 into master 
integrals. But if we also consider the previous identities of the group 3, we find one relation which 
transforms the whole combination W2- The presence of this relation makes it unnecessary to search for 
identities which singly transform the integrals of W2. 



3 Difference equations for master integrals 

By using the algorithm shown in the above sections, a generic Feynman integral can be reduced to master 
integrals. Now we consider their calculation. 

We consider the generic integral (7) as a function of the exponents U{a,(})] the integration-by-parts 
identities (6) form a system of recurrence relations satisfied by the function U. The key observation is 
that such recurrence relations are linear difference equations^ satisfied by the function U{a,P). This 
observation has two important implications: 

1. Theory of linear difference equations[10] is a well established (but not very well-known) mathemat- 
ical topic for which various useful mathematical tools exist. 

2. Difference equations can be solved numerically, and therefore this establishes a new method of 
calculation of the integrals U{a,(3). 

Another important observation is that the integration-by-parts provides recurrence relations in several 
variables for the functions U, that is, a system of partial difference equations. In general, the numerical 
solution of such a system is, from the point of view of the numerical calculation, a task comparable to 
the solution of a system of partial derivative equations or to the multidimensional integration, so that 
high-precision results may be very difficult to obtain if the number of variables is high. Therefore, if we 
desire to obtain high-precision results, we must consider only difference equations in only one variable. 

In section 3.1 we show where to introduce such single parameter; in section 3.2 wc show how to 
obtain recurrence relations (that is, difference equations) in this parameter from the integration-by-parts 
identities (6). 

3.1 Difference equations in one exponent 

Let us consider, for instance, a master integral of the form 

^-J D,D2...D,^ ■ ^''^ 
We raise one denominator to power x. If we choose Di, we define a function Udi 

UnAx) = J DID2...D,, ' 

where the subscript shows that the exponent has been introduced in the denominator D\. The value of 

the integral (18) is recovered as Udi (1) = B. 



® Other authors are also encountering other kinds of difference equations in the evaluation of diagrams [11, 12] and 
solving them, by using techniques different from the standard methods of [10] used in this paper. 
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Using intcgration-by-parts identities, we look for an identity which relates Udi with integrals with a 
smaller number of denominators. As i? is a master integral, an identity which directly transforms Udi {x) 
into simpler integrals cannot be found; instead, as we will see in the next section, one finds identities of 
the form 

R 

Y,Vj{^)UDA^+3)=F{x) , (20) 

3=0 

where pj {x) are polynomials in x, and F{x) is a known function. U Oi appears with arguments shifted 
by an integer. This identity is a linear difference equation of order R in the variable x satisfied by the 
function Udi{x). 

The right-hand side function F{x) is in general a linear combination of functions analogous to Udi , 

derived from master integrals containing Di, but with some of the denominators D2. -D3, . . . , missing. 

We note that it is possible to raise to x another denominator, D2, D3, ... , etc. (or equivalently to 
permute the numbering of the lines of the diagram); for each denominator Dj a function Uojix) will be 
defined. In general the functions U jj. (.x) arc different each one and satisfy different difference equations, 
but for X = 1 they all have the same value ?7_d^. (1) = B. This fact is particularly useful for checking the 
consistency of the calculations. 



3.2 Construction of systems of difference equations 

Let us suppose that, by means of the procedure described in section 2.5, we establish that there are L 
master integrals, with a number of denominators ranging from Nk to A^^- 

We note that some master integrals may have the first denominators missing. Then wc; group the 
master integrals in Af^ — TV^ + 1 sets Sm such that each set includes the integrals which contain the 
denominators Di with i> m. Each set will contain master integrals Bmi, 

[rf^fcl]...[d^fcjvj ^ , m<i2<...<in, (21) 

where m = 1, . . . , A^^ — A^^; + 1, / = 1, . . . , and = L; the number of denominators n ranges 

from Nk to Nd — m + 1. Note that n, the indices i and the exponents /3 of the various integrals depend 
on m and /. The ordering B„ii, Bm2, ■ ■ ■ , -SmL„ follows the ordering of master integrals with the same 
first denominator Dm provided by the algorithm 1. Changing the exponent of Dm from 1 to x in each 
master integral, we define the "master functions" 

[d^fci] . . . [d^fc^vj ; (22) 

Umi{^) = Bmi recovers the original master integral. Now we must find the difference equations satisfied 
by the functions Umi{x). We follow a method quite similar to that used for the reduction of a generic 
integral to master integrals. We build a system of identities obtained from integration-by-parts using 
Eq.(6), modifying the integrand V^^^p of Eq.(5) by adding x — 1 to the exponent of the first denominator: 

In this way all the integrals become functions of x 

Unia0{x) = J [rf^fcl] . . . [d'^kM,] Ki„^(x) , (24) 

and the identities (6) become difference equations between these functions. Integrals which have a 
different first denominator D^ never appear in the same identity, so that it is convenient to build and 
solve separately the systems of the identities with different values of m. We do this with the algorithm 

Algorithm 2 Consider the algorithm 1 with the following modifications: 

1. Set i\ = m. 



10 



2. Add X — 1 to the exponent of Dm ■ 

3. Ignore the addition of x — 1 to M4 in the definition of set (9). 

4- Add the following conditions to the order of priority of the extraction of integrals: 

(a) the generic functions Unia0 must be extracted before of the functions Umi; 

(b) the functions Umi{x — 1) (generated by the algebraic identities) must be extracted before of the 
functions Umi{x + i), where i > 0. 

As result of the solution of the system, we find lots of identities analogous to Eq.(lO), not interesting, 
which express generic functions Uniap as combinations of master functions Umi with argument possibly 
shifted, and rational functions of x, D, masses and scalar products of external momenta as coefficients; 
we also find much more interesting identities containing only master functions Umi- 

Let us suppose, for simplicity, that each master integral Bmi contains a different combination of 
denominators. With a suitable choice of the constants and bi, the solution of the system provides 
difference equations satisfied by the functions Umi, with the structure 

Ri l-i Qik 

^Pii{x)Umi{x + i) = ^^qjki{x)Umk{x + j) , 1 = 1,... ,Lm , (25) 

j=0 fc=l j=0 

where pu and qjki are polynomials in x, D, masses and scalar products of external momenta. The 
right-hand side of the ^th equation contains the functions from Umi to Um,i-i- Therefore the set of L,,, 
equations (25) forms a triangular system of difference equations. The triangular structure is particularly 
useful for simplifying the numerical solution: the equations are solved in ascending order, one at a time, 
for Z = 1, 2, . . . , Lm] when the equation for Umi has to be solved, the functions Umi, ■ ■ ■ ,Um,i-i in the 
right-hand side are already known. Note that some functions may be missing from the right-hand side; 
for example, if the master function Umi has Nk denominators, the right-hand side is zero. 

Now wc consider the more general case where G different master integrals have the same denominator 
{Di^---Di^) and different numerators. These master integrals correspond to master functions with 
contiguous indices Um,i'+i, Um,i'+2, ■ ■ ■ , Um,i'+G as effect of the ordering of master integrals. The 
solution of the system provides for these functions a subsystem of simultaneous G difference equations of 
the kind 

G Rhs r Qhi'k 

Y.Y.P'wl'h{^)Um,l'+g{x + l)=Y.Y.l'okl'h{x)Urak{x+j), /l=l,...,G, (26) 
g=l i=0 k=l 3=0 

not in triangular form, where the left-hand side of each equation contains all the master functions Um,i'+\, 
. . . , Um,i'+G- We prefer not deal with the solution of the subsystem of simultaneous difference equations, 
so that we transform it into triangular form; this is not obligatory, is only convenient to simplify the 
subsequent numerical solution. We make the replacement x ^ x + c in the equations (26), with c = 
1,2,..., generating new identities which are inserted in the subsystem, and taking care of extracting the 
terms containing the function Um,i'+j before of the terms containing Um,i'+k if j > k. This procedure is 
repeated until one obtains a set of G equations in triangular form, 

K I'+h-iQ'hi'k 

J2Pil'hix)Um,l'+h{x + i)= J2 J2 Ijkl'hixWmkix + j) , h=l,...,G. (27) 

1=0 k=l 3=0 

Another advantage of the transformation into triangular form is that one obtains a difference equation 
containing the function Um.i'+i, but not containing the other functions Um,i'+2, Um,i'+3, • . • , so that 
Um.i'+i, which corresponds to the master integral with 1 as numerator, can be found independently. 
Unfortunately the transformation into triangular form has a price: one obtains for Um.v+i an equation of 
order R'^ > maxh Rhi which has coefficients q'Jkii lix) much more complicated and cumbersome than the 
coefficients qj^fhix) of the equations of the subsystem (26). If R[ is large (say R'l > 10) these coefficients 
become huge and difficult to obtain so that it may be more convenient in these cases to solve the system 
of simultaneous equations (26) directly; this will be described in some next paper. 



11 



Once all the subsystems of simultaneous equations corresponding to the various groups of master 
integrals with equal denominators are transformed into triangular form, the whole system takes the form 
(25). 

Concerning the choice of a, and 6,, all the considerations done in section 2.6 remain valid; for each 
combination of n denominators {-Di^ . . . Di^} there is a minimal subsystem n: Q "g' whose solution 

allows one to obtain the equations (26). For the diagrams so far examined di turns out to be always 
equal to the maximum number of scalar products appearing in the master integrals (typically 1 . . .3); 
typical values of bi are 0, 1, 2. 

4 Solutions of difference equations with factorial series 

By using the algorithms of section 3, the triangular system of difference equations satisfied by the master 
functions is worked out. Now we consider the solution of each difference equation. 
Let us suppose, for instance, that the master function U{x) defined as 

U{x) = J [d^h] . . . [d^fc^vj ^^x^, j^^^ ; (28) 
satisfies a difference equation of order R 

R 

Y,Pi{x)U{x + i) = F{x) , (29) 

where Pi{x) are polynomials and F{x) is some known function. The solution of this nonhomogeneous 
equation can be written as 

U{x) = U"°{x) + U^"{x) , (30) 

where U^^ is a particular solution of the nonhomogeneous equation (29) and is the general solution 
of the homogeneous equation 

R 

Y,Pi{x)U"°{x + i)=0 . (31) 
The general solution of Eq.(31) can be written as 

U"0{x)=Y.^,{x)Uf''{x), (32) 

where LOj{x) arc periodic fimctions of period 1, and {U^'^{x), .... U^'^{x)} is a fundamental system of 
independent solutions of the homogeneous equation. In the following, recalling from [10] for the ease 
of the unfamiliar reader the essential matter on the solution of linear difference equations with factorial 
series, we describe how to obtain factorial series expansions of C/^*^ and U^^. 

4.1 Factorial series 

The series 

f^^T{x-K + s + l) T{x- K+l)y'°^ X- {x - K + l){x - K + 2) 7(33) 

is known as factorial series of the first kind [13], or series of inverse factorials [14], or faculty scries [15]. 
We refer to it in brief as factorial series. The series converges for every point in a half-plane which is 
limited on the left by iRx = A (excluding x = K — 1, K — 2, . . .). The number A is the abscissa of 
convergence. If A = 00 the series is everywhere divergent. 

As coefficients of the series encountered in this work behave for large s as |as| ~ s!s", it is useful 
(especially for numerical applications) to define the reduced coefficients a'^ ~ as/s\. For large s the 
generic term of the sum tends to a'gT{x + l)s^~^, so that the factorial series has the same convergence 

Lo 



properties as the Dirichlet series Yl'^o ^'s^^ ^ 
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4.2 Operators tt and p 

Given an arbitrary number m one defines the operator p as^ 

p-Ui.) = ^^^^±^^Ui.-m). (34) 

This operator has the property 

p"'p''U{x) = p™+"f/(a;) ; (35) 
if the operand is the unity, we omit it and write 

p"»l = p"»= ^(^ + ^) ■ (36) 
r(a; — m + 1) 

One defines the operator tt as 

TTUix) = x{U{x) - U{x - 1)) . (37) 
The following properties can be proven: 

[K,p\U{x)=pU{x) , 

p{tt)p^U{x) = p'^pin + m)U{x) , ^ ' 

xU{x) = (tt + p)U{x), 

x-U{x)=f2[f2{-iy-^(As,ui^--^\ p'^Uix) , (39) 

p{x)U{x) = [p(7r)+pi(7r)p + p2(7r)pV2! + .--] U{x) , 
where p is a polynomial, Pn is 

p„(A) = A„f,(A) = f^(-l)^ (fjpiX - J) , 

and Sjk are the Stirling's numbers of second kind [16]. Using these operators an expansion in factorial 
series becomes an expansion in powers of p~^ so that we will able to obtain solutions in factorial scries of 
difference equations in the same manner as power series solutions of differential equations are obtained. 

4.3 Solution of the homogeneous difference equation 

Let us consider the homogeneous difference equation (31) of order R 

Po{x)U"'^{x)+pi{x)U"'^{x + 1) + ...+pr{x)U"''{x + R) = 0; (40) 

making the replacement x ^ x — R and defining qii-i{x) = Pi{x — R) the equation becomes 

qo{x)U''^{x) + qi{x)U"0{x - 1) + . . . + te(a;)[/^0(a; -R) = 0. (41) 

Using the previous definitions of operators now wc search for a formal solution in factorial series. We 
make the change of variable U^^{x) = ji^V^^ix) in Eq.(41), where n is an unspecified parameter, 
obtaining 

li^qn{x)V"°{x) + n"-^qi{x)V"^{x - 1) + . . . + qR{x)V"° {x - i?) = . (42) 



In Ref.[10] the operators tt and p are defined more generically as p™U{x) = (r{x — r + l)/r{x — r — m + l))U{x — m) 
and TrU{x) = (x — r){U{x) — U{x — 1)), where r is a fixed number. Here for simpUcity we have set r = 0. 
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Now we multiply the equation by x{x — l){x — 2) . . . [x — R + 1) and we observe that xV{x — 1) = pV{x), 
x{x — l)V{x — 2) = p'^V{x), etc.; the equation takes the form 

[(/.o(a;,M) + <^i(x,/x)p + .-. + <^fl(a;,M)p'*] V"^{x)=0, (43) 

where (f>j arc polynomials in x and fi. The miiltiplication by x is equivalent to the multiplication by 7r + p, 
therefore substituting the relations (39) in Eq.(43) one obtains the first canonical form of the difference 
equation: 

[fo{7v, fi) + /i(7r, fi)p + hin, fi)p' + ... + /^+i(7r, m)p'"+'] V'^'^ix) = , (44) 

where fi are polynomials in tt and /x. In the case of the difference equations encountered in this work 
we find that fm+i{''^, jj) is independent of tt, but not independent of fjb; therefore /m+i(7r, jS) = fra+i{lj)- 
The algebraic equation in 

/„+i(/x)=0 (45) 

is the characteristic equation^. It turns out that the characteristic equation has always R solutions 
different from zero. Let /ii, /i2, • • • , Ma be the A distinct solutions of Eq.(45). For each of these values, 
/X = yUj, i = 1, . . . , A, the first canonical form (44) takes the form 

[/o(7r) + /i(7r)p + /2(7r)p^ + . . . + /™(7r)p'"] V^'^ix) = . (46) 

Now we try to satisfy the equation (46) in V^*^ with the factorial series 



\rHO(^\ _ \^ ggi \X + ij \ ^ , K-s _ „ _u„ 



^T{x-K + s + l) 

s—0 s—0 



+ ..., (47) 



whose asymptotic behaviour for large x is V^^{x) ~ aox^ . Substituting Eq.(47) in Eq.(46) one obtains 

the recurrence relations 

ao.frn(K + m) = , 

aifrniK + m - 1) + aofm-i{K + m - 1) = , 

asfm{K + m- s)+as-ifm-i{K + m- s) + ... + as-mfo{K + m- s) = (s > m) . ^^g^ 

Supposing that oq 7^ 0, the equation 

fmiK + m)=0 (49) 

is the indicial equation. Let Ki, K2, ■ ■ ■ , Ki, the roots of this equation. In the case of the difference 
equations encountered in this work all the roots turn out to be distinct, and v turns out to be the 
multiplicity of /i,. For each of these values of K we solve the system of recurrence relations. If there 
arc no roots differing by a positive integer, ,fm{K + m — s) 7^ for s > 1, and therefore the recurrence 
relation (48) provides the coefficient as of the factorial series for every s. If there are roots differing by 
a positive integer (so-called congruent roots) we find fm{K + m — so) = Qiov some sq, so that the term 
0'sa,fm{K + m — .So) vanishes from the relation (48) used to obtain as„. The remaining terms of this 
relation always vanish^ in the case of the difference equations encountered in this work, therefore the 
value of remains undetermined and can be chosen at will (usually one puts = 0). For details on 
the convergence of the factorial scries so obtained, sec section 6.1. 

In order to obtain the general solution of the difference equation, for each distinct solution of the 
characteristic equation /Xj, i = 1,... ,A, we find an indicial equation, whose solutions are Kij, j = 

* If the coefficients of the equation are Pi(x) = '^^j—QPij^^~'' t the characteristic equation has the explicit expression 

Otherwise V^'-^(x) would be a derivative with respect to if of a factorial series, with asymptotic behaviour 
~a;^log"(i;). 
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1, . . . , Mi. For each solution of the indicial equation we solve the system of recurrence relations between 
the coefficients ai*'''^ and we find Vi solutions of Eq.(46), 

Multiplying by ji^ we find solutions of the difference equation (40). In total we find Yl^=i Ui = R different 

solutions. The general solution of the homogeneous difference equation (40) will be a linear combination 
of the these solutions with periodic functions tUij{x) of period one as coefficients, 

U"0{x) =Y.Y.Q,,{x)„^V,^°{x) . (51) 
4.4 Solution of the nonhomogeneous difference equation 

Let us consider the nonhomogeneous difference equation (29) of order R 

p^{x)U''"{x)+p^{x)U''"{x + l) + ...+pn{x)U''"{x + R)=F{x) ; (52) 
making the same replacements as Eq.(40) the equation becomes 

qo{x)U''"{x) + q^{x)U''"{x - 1) + . . . + qR{x)U''"{x - R) = F'{x) , (53) 

where F'{x) = F{x — R). In the difference equations encountered in this work, the right-hand side can 
be written as a sum of factorial scries expansions. We assume for simplicity that F'{x) = ii^T(x) and 
T{x) has the factorial series expansion 

T{x) = CqP^ + cip^-i + C2P^-2 + . . . , (54) 

where K and Ci are known. 

Making the substitution U^^{x) = ^^V^^{x) in Eq.(53) and using the relations (38-39) we obtain 
the canonical form 

[/o(7r) + /i(7r)p+ /2(7r)p2 + . . . + /„(7r)p'"] V^'^ix) = T{x) . (55) 

The expansion in factorial series of V^^{x) has the form 

V^'^ix) = aop""-"" + aip^-™-i + aap^-"-' + . . . ; (56) 

inserting this expansion in Eq.(55) and equating the coefficients of the powers of p one finds the recurrence 
relations 

a.ofm{K) = Co , 

aifm{K + aofm-i{K - 1) = ci , 

(57) 

asfrn{K - s)+ as-ifm-l{K - s) + . ..+ as-m.fo{K - s) = Cs (s > m) . 

Solving this system we can find the coefficients a,, . If /,„ {K — sq) = for some .so ■ the term /,„ (K — sq) 
vanishes from the relation (57) and, in the difference equations encountered in this work, the remaining 
terms form the trivial identity Cg = Cg, so that the value of a^g remains undetermined and can be 
chosen at will. If the factorial series (56) converges, U^^{x) = ji^V^^ {x) is a particular solution of the 
nonhomogeneous difference equation. 

The case F'{x) = n^p{x)T{x), where p is a polynomial, can be reduced to the previous case by 
transforming x into tt -|- p in the polynomial and letting the operators act on the expansion of T{x) 
following Eq.(38). 
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5 Determination of arbitrary constants 



The general solution of the difference equation (29) can be written as 

R 

U{x) = ^w,(x)i7f °(x) + U^"{x) , (58) 
i=i 

where CiJj{x) are periodic functions of period 1. If we consider only integer values of x, the values of(jjj{x) 
are independent of x so that we can replace them with arbitrary constants r]j : 

R 

U{x) = Y^rjjUf'^ix) + U^"{x) . (59) 

The unambiguous solution of the difference equation requires the determination of the R constants 
rjj. The value of these constants can be determined 

1. from the large-a; behaviour of the solution (59) and of the integral (28), by equating the first 
coefficients of the expansions in factorial series; 

2. by equating the values at a; = of Eq.(59) and Eq.(28). 

As wc will sec in the following sections, the method 1 may provide the values of all the constants, 
but the determination of the coefficients of the factorial series turns out to be a simple task only for 
euclidean massive integrals (see sections 5.1, 5.2 and 5.3), where it involves integrals with one-loop less, 
and gets more complicated for non-ouclidcan integrals or integrals with zero masses (sec section 5.4); on 
the contrary the method 2 has the limitation that it provides only one relation between the constants, 
but on the other hand it is not affected by the kind of the external momenta or masses (see section 5.5). 



5.1 Lcirge-a; behaviour of integrals: euclidean massive case 

In this section we work out a relation between the coefficients of the expansion in factorial scries of the 
master integral (28) and simpler integrals with one-loop less. We assume that all the external momenta 
are euclidean, that is, the matrix of the scalar products pi ■ pj is semidefinite non-negative, and no mass 
is zero. 

Choosing the momentum routing of the integral (28) such that Di = kl + mf, we write it as 

where g includes the contribution of the remaining Nk — 1 loops 

gik,) = [ [d-k,] . . . [d-k,,] "^"^r^^^'";!"'"-^^'^ • (61) 
J 1J21J3 ■ ■ ■ 

Note that the function g also depends on the external momenta pj and the masses m2, ma, ... , mjv^ of 
the other denominators. From a graphical point of view, g{k\) corresponds to the original diagram with 
the line Di cut. Introducing hyperspherical polar coordinates for the integration over ki 

d'^ki = \kif-'^d\ki\ dflnik) , (62) 

[Qd = 27r^/^/r(D/2) is the £>-dimensional solid angle) and separating the angular and radial part, 
Eq.(60) becomes 

00 

1 r rlh"^ (U2\D 2-1 



where / is the angular mean over ki of g 

f{kl) = -^j dilnik) g{ki) . (64) 
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For large x the factor {k1+m\)~^ of Eq.(63) peaks strongly around k\ — 0; because of the assumption 
on external momenta and masses, the function /(fcj) has no singularities for ki > 0, and therefore we 
expect that the large- a; behaviour of the integral is controlled only by the behaviour of f{k\) near k\ = Q. 

Making the change of variable ki = rui- in Eq.(63) we obtain 





where f{u) = f{m\u/{l — u)). We expand f{u) in u 



~f{u) = u''{l-ufY.^su'; (66) 

s=0 

a is an integer greater than or equal to zero (/ is regular in the origin), while the factor (1 — uY has been 
introduced for convenience. Expressing the integrals over u in terms of Beta function (which takes care 
of analytical continuation if the integral diverges at the endpoints), and choosing (3 = D /2 + 1, Eq.(65) 
takes the desired form of a factorial series 

^(.) = ,g^a.p^o-.^,g^,^ (^_+^^ , (67) 

s=0 s=0 ^ ' 

where 

lio = l/ml, Ko = -D/2-a, and = b^mf r(s + £>/2 + a)/r(£)/2) . (68) 

The coefficients bg of Eq.(66) can be easily expressed in terms of the coefficients of the expansion in 
kjoifikl) 

OO 

/(fcD = (fclrE/^^l^ (69) 

{bo = mf"/o , ■ ■ ■), so that the largc-x behaviour of U{x) proves to be determined by the behaviour of 
f{ki) for small kf. For large x the leading behaviour is given by the first term of Eq.(67) 

In the frequent case of an integral with numerator equal to one, a = and /o = /(O), it becomes 

U{x) « (m?)^/2-^a;-^/V(0) • (71) 

Now we equate Eq.(67) with the expansion in factorial series of the general solution (59) (see Eq.(51) 
and Eq.(56)), obtaining 

In this equation only the R constants r]j are unknowns, each one corresponding to a solution Uj^*^ with a 
different pair of values fXj and Kj , fij being one solution of the characteristic equation of the homogeneous 
equation and Kj being the solution of the corresponding indicial equation; we have set djo = 1. We have 
also assumed that the expansion of the nonhomogeneous term contains a sum of expansions with different 
pairs of values ^if^^ and Kf^^ . 

Fortunately, the number of constants rjj to find can be drastically reduced. If the difference equa- 
tion is homogeneous, the constants rjj which may be different from zero are only those such that the 
corresponding /Xj and Kj satisfy the condition 

lij = l/m\, Kj+D/2 + a = mteger<Q ?7j 7^ ; (73) 
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all the other constants must be zero. 

In the case of nonhomogeneous equation, one must recall that it is part of a triangular system of 
difference equations; the nonhomogeneous term receives contributions from master integrals with a smaller 
number of denominators, which in their turn satisfy other nonhomogeneous difference equations, up to the 
simplest master integrals which satisfy homogeneous equations. Clearly jj,^^ = l/m\ and Kj^^ + D/2 
is always an integer, therefore rjj may be different from zero if fij and Kj satisfy Eq.(73), or the new 
condition 

lij = l/m\, Q <Kj+D/2 + a = integer < mscKKj^" + D/2 + a % 7^ (74) 

which implies cancellations between the first coefficients of the homogeneous and nonhomogeneous ex- 
pansions. We found that in the equations considered in this work, the condition (74) is never satisfied, 
and the condition (73) is satisfied by a small number of pairs of /ij and Kj, often only one; if no pair 
satisfies the conditions, U {x) is completely determined by the nonhomogeneous term. 

The values of the non-zero constants Tjj are determined by comparing the coefficients of the same 
powers of p of the two sides of Eq.(72). The required few coefficients as are easily calculated from the 
first coefficients fs of the expansion in k\ of f{k\). 

The required coefficients /« are calculated by expanding numerators and denominators of Eq.(61) for 
small ki , and by performing the angular integration over ki of Eq. (64) . Angular integrals are straight- 
forward as they contain exclusively powers of scalar products containing fci. As a result, the coefficients 
fs will be expressed by integrals over fc2, . . . kfq^, belonging to diagrams with one loop less (for example, 
if the numerator of Eq.(61) is unity then /o = /(O) — .9(0)). These integrals can be expressed using 
algorithm 1 as combinations of the master integrals of the new diagrams; in their turn, these new master 
integrals can be calculated by inserting an exponent in a denominator and building and solving new 
difference equations. The arbitrary constants of their solutions can be found using the large-exponent 
behaviours, which can be expressed in terms of integrals with another loop less, and so on. In this way 
we can explicitly calculate the values of all the arbitrary constants r/j. This fact is very important. 

5.2 Non-euclidean case: deformation of the radial path 

Now we consider the non-euclidean case: the matrix of the scalar products Pi ■ pj is definite negative. We 
write Eq.(63) as 

lo 

where we explicitly show the dependence on the external momenta P = {pi, . . . ,PNp}, and where Iq 
indicates the path of the radial integration. For the moment the path is assumed to be the positive axis 
in the complex plane. 

The integral U{x,P), considered as a function of the external momenta P, is defined in the non- 
euclidean region by an analytical continuation through a generic path in the complex P-space. The path 
begins in some euclidean initial point Pin and ends in the desired non-euclidean final point Pend- It is 
possible that in some intermediate point of the continuation path in the P-spacc, for a particular value of 
the external momenta P,, /{kf, Pg) is singular in a point k\ = > 0, just on the radial integration path, 
breaking off the analytical continuation. If this singularity in ki cannot be avoided by modifying the 
continuation path in the P-spacc, then we are forced to deform'^'^ the integration path in the fc^-spacc by 
turning around the singularity; the singularity in k\ itself moves in the complex /c^-plane to a final point 
k1^^ when we complete the analytical continuation from P, to Pend- The final integration path Iq starts 
in ki = 0, turns around the singularity fc^^^^ (usually on the negative axis) and comes back to kf +00. 
In the general case of multiple singularities the path turns around all them. The deformation of the 
radial path must be performed when the values of the scalar products pi ■ pj get over some "deformation 
thresholds" in the non-euclidean region, the exact point of these thresholds being depending on the values 
of masses and the structure of the diagram. 

One can show using Feynman parameters that the deformation thresholds are the thresholds (anoma- 
lous thresholds included) in P of g{0,P), the diagram obtained by eliminating the line Di and setting 

1" The deformation of the radial integration path was first discussed in [17] for self-energy diagrams in D = 4. 
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ki =0 everywhere. Possible denominators depending only on ki give "additional" deformation thresholds 
of the kind = — m^. 

Once the deformation thresholds are determined, we can split the non-euclidean region of the P-space 
into two regions: 

1. one region below the deformation thresholds which adjoins the euclidean region, where the integra- 
tion path is the same as the euclidean region; 

2. the remaining region above the deformation thresholds where the integration path turns around 
some singularities. 

5.2.1 Example: deformation of the path for the one- loop self-energy integral 

As an example, now we consider the simplest case, the one-loop integral 



J {e+ml){{p-kf+ml)~ T{D/2)] k'^+ml nnJip- 



kf + ml 



(76) 



The angular integral is 



11 



1 f dnnik) Z^^z')^ (77) 



J {p — kY+m\ pk 
where ^(Z^) is the hypergeometric function F{1, 2 - D/2; D/2; Z"^), 

^^,^ p^ + e+mlTR{p\k\-ml) ^ ^^^^ ^^^^ ^^^^ 

and R{x, y, z) is the usual two-body phase space square root 



R{x, y, z) = ^x^ +y'^ + z'^- 2xy - 2xz - 2yz . (79) 

The angular integral (77) has two branching points in the complex fc^-plane, fc^ = {p±im2)^. Considering 
Eq.(76) in the euclidean case, > 0. the path of integration is always the positive real axis. If < 0, 
Eq.(76) must be analytically continued in the complex p— plane from a generic initial euclidean point 
P = Pin > to the non-euclidean final point p = i\/ —p^ which is on the imaginary axis. The initial point 
has 9pi„ = 0. If p^ > — m2 no radial singularity appears so that 

^-fiDj2)J k^+ml pk^^^>- 



If p^ < —rn^, the final point has 3p > TO2 so that any path connecting the initial and final point must cross 
the line 3p = m2 in some point. Let pc = im2 + c be this point; if we consider Eq.(77) with momentum 
Pc, one of the branching points falls vak^ — > which is exactly on the path of integration. Therefore 
the path of integration must be deformed avoiding the branching point (p — 11712)^ which is on the 
negative real axis. One chooses as path of integration [17] first the segment ((p — zto2)'^,0) above the 
negative real axis, taking the square root R with the plus sign and using in the place of Z, then the 
straight line ((p — im2)^,oo) below the axis. Therefore, if p^ < —m\ Eq.(76) becomes 

(p-im2)^ , 00 

dfc2(A;2)^/2-l Z-1 1 r dfc2(fc2)^/2-l Z 



J k-^+ml pk ^ ' T{D/2) J k'^+ml pk ^ ' ' 



r(Z)/2) J i\, r iii'i pry i.yjj/z,j j n,- -r nii yn, (Q^\ 

This result will be used in section 7.2. 



Here we generalize to arbitrary D dimensions the result obtained in [17] in the four-dimensional case. 
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5.3 Large-x behaviour of integrals: non-euclidean massive case below cind at 
the deformation threshold 



In the region below the deformation thresholds, in the case of non-euclidean massive integrals, the deter- 
mination of the large-x behaviour is quite similar to that described in section 5.1. 

At the deformation threshold the determination of the large- .t behaviour presents some peculiarities, 
since some external momenta have on-mass-shell values. In the integral (61) denominators of the form 
kf — 2p ■ ki appear, where p is some on-mass-shell momentum. Considering the expansion in fci of g{ki), 
these denominators vanish at ki = and cannot be expanded in a way as straightforward as for off- 
mass-shell denominators; they must be included in the angular integrals, complicating quite a lot the 
integration. 

In order to explain the situation, let us consider here the case of a one-loop diagram with N + 1 
external lines, and the integral 

[d'^k] 



(fc" + m^)((pi - ky + mi)... ((pat - ky + m%) 



(82) 



Setting = we see that the deformation threshold is pf = — m?, i = 1, . . . , A'^. Now we set the external 
momenta to these on-mass-shell values, and we consider the integral 

/ [d^ k] 

^^^^ " J (fc2 + m2)^(fc2 -2pi-k)--- (fc2 - 2pN ■ k) ' ^^^^ 

whose we want to calculate the large-x leading behaviour. Following the notation of section 5.1 we write 



W{a 



_ f dk\er/--^ m 



f(k') = ^ / ^^^(fe) (85) 

no J {k^ -2pi-k)---{k^ -2pN -k) ' ^ ' 

Now we extract the leading behaviour of /(fc^) for fc^ — > 0. Inserting Feynman parameters Xi, i = 1, 
. . . , N — 1 one finds 

T(N) f f dnoik) , , 

where P{xi) = J2iLi ^iPi ^nd J2iLi = 1- The angular integral in D dimensions can be expressed using 

a generalization of the formula (77) 

' I' ^"^^^^ ^^yFix,x + l-D/2-D/2-Z'); (87) 



J {{p — ky + m?y \pk J 

using properties of hypergeometric function one finds for k'^ ^ 

J_ f dnnjk) ^ , r'(x)h')-^^' TiD/2)T{N/2) 

no J (fc2 - 2Pixi) -k)^ ^ ' 2V{N)T{{D - N)/2) ' ^ ' 



then 

^ nDimN/^i, 2.-N/2 f dx^...dxN., 

^^"^^ 2T{{D-N)/2f > J (-P2(^^^iv/2 • 



We see that f{k^) is singular in fc = and, if N > D, the integral over k^ of Eq.(84) diverges for small 
k^; Eq.(70) with arbitrary a takes care of the necessary analytical continuation. 

Let us now consider a one-loop integral in E dimensions with the same denominators as the angular 

integral (85) 



[d^q] 



9) • • • (g2 _ 2pN ■ q) 



(90) 
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introducing Feynman parameters in the same way as for Eq.(86) we obtain 



If we choose E = N the integrals over Xi of Eq.(89) and Eq.(91) become identical, so that we can rewrite 
Eq.(89) as 

f('=')-('^r'" \^^^'^%„^ L„M. (92) 

This result shows that the angular integral (85) with N on-mass-shcll denominators in the limit ^ 
is proportional to a one- loop integral with the same denominators (one less than Eq.(83)), calculated in a 
number of dimensions equal to the number of on-mass-shell denominators. Moreover, using Eq.(92) and 
Eq.(70) (with a = —N/2 according to Eq.(89)) one finds the large-x leading behaviour of Eq.(83): 

W{x) « (m2)^/2-^/2-a;-^/2+^/2iw(p,) . (93) 

This result is similar to euclidean leading behaviour for a one- loop vacuum integral (Eq.(71) with /(O) = 

1); the differences are the change of the exponents and the multiplication by the factor Ljv(k) which 
is independent of the number of dimensions D of the integral (83). As a consequence of the divergence 
of the radial integral and of the necessary analytical continuation, we observe an apparently paradoxical 
result: if A'' > D and uiq = 1, from Eq.(93) we see that the integral W{x) increases as the exponent 
X increases, while in the euclidean case the integral U{x) decreases as x increases. The integral L^ipi) 
may be calculated analytically, for example, extracting the A;^ — » behaviour from the known analytical 
expressions of the angular integral (85) for = 4 and A'' = 1 [17] and = 2, 3 [18], or numerically, by 
inserting an exponent in one denominator of Eq.(90) and building and solving a difference equation, or by 
using the identities of the section 5.5 (for example see Eq.(152)). In section 9.2 we will need the following 
values Ln of L^lpi) in the on-mass-shell and equal masses case m, = 1, = — 1 and {pi — Pj)^ = — 1 
for every i and j: 

Li =V^/2 , L2 = ttV3/9 , 

/ r- r-\ , ~ (94) 

L-i = ( arctan V8 - arctan VS j V97r/8 , L4 = 0.172751462 .... 

The generalization of Eq.(93) to multi-loop integrals is straightforward; considering 

/[d^k] 



the large- a; leading behaviour is 

W'{x) « (m2)^/2-^/2--a;-^/2+iv/2j^^(p.^ ^(0^ _ (97) 

The condition (73) in this case must be modified to 

= l/mf , Kj+D/2- N/2 = integer or half-integer < VjT^O ■ (98) 



5.4 Non-euclidean case above the deformation threshold and massless case 

The cases described above, euclidean and non-euclidean below the deformation thresholds, are both united 
by the fact that the large- a; behaviour of U {x) depends only on the behaviour of /(fci) near the origin. In 
these cases it is easy to express this behaviour in terms of simpler diagrams obtained by putting fci = in 
the main diagram. The situation gets more complicated if f{k\) has other singularities on the integration 
path: for the ith. singularity placed in g? 7^ 0, an additional contribution to the large- a; behaviour of U{x) 
appears, of the form (g? -|- m\)~'^V{x), with V{x) x^ for large x. The precise form of V{x) depends 
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on the behaviour of / near the singularity; unfortunately its determination requires the calculation of 
angular integrals for kf^O (see Eq.(64)) 

J driDik) g{ki) , (99) 

which is difficult and case-dependent. In this work we have performed these calculations only in the 
one-loop example of section 7.2. In the case of more complicated diagrams we have preferred to avoid 
the calculation of these angular integrals and to look at the problem from a different point of view, using 
a method with broader applicability based on differential equations (see section 10). 
The cases where singularities of / other than the origin appear are: 

• The non-euclidean massive case, above the deformation threshold: the path of integration must be 
deformed, turning around a number of singularities of /{kf). See section 7.2.1 for an example. 

• The zero mass case, where mi > and some of the masscis m2, . . . , rriN^ arc zero. If the external 
momenta are euclidean, singularities on the positive axis of kf may appear, but this fact does not 
require deformation of the path. See section 7.2.2 for an example. 

The situation worsens if the mass mi of the denominator raised to x is zero. In this case the behaviour 

of U{x) for large x depends on the behaviour of f{kf) on the whole real axis, and not in some isolated 
points. This can be easily understood by noting that in this case U{x + D/2 — 1) is the Mellin transform 
of f{kf). See section 7.2.3 for an example. 



5.5 Zero exponent condition 

A very useful relation between the constants rjj arises from Eq.(59) calculated at x = 0: 

R 

U{0) - U^"{0) = ^r;,i7f 0(0) . (100) 

?7(0) is an integral without the denominator Di] if the value of U{0) can be found by some method, as an 
example, by adding an exponent to D2 and by solving a difference equation, and if at least one Uj^'-'{0) 
has a non-zero value, Eq.(lOO) establishes one new relation between the constants r]j. The existence of 
this relation must be verified in each particular case. The advantage of this relation over the relations 
found using the asymptotic behaviour of U{x), is that it is valid for every kind of integral, regardless of 
the value of masses or external momenta. 

It is possible to construct identities analogous to Eq.(lOO), by inserting in Eq.(59) x = —1, —2, . . . 
instead of a; = 0. Unfortunately in all the analyzed cases the identities so found turned out to be 
equivalent to Eq.(lOO). 



6 Evaluation of factorial series 

Once the constants r]j of Eq.(59) have been determined using the methods described in section 5, in order 
to obtain the value of the master integral U{1) we must calculate the values of the homogeneous solutions 
Uj^^{l) and of the nonhomogeneous solution U^^{1). 

6.1 Convergence of factorial series and instabilities of recurrence relations 

Let us unify the notation by considering the solution U^"\x), where (a) indicates one of the solutions of 
the homogeneous or nonhomogeneous equation, and expand it in factorial series: 

Convergence of the scries depends on the value of the abscissa of convergence A. Analyzing the large-s 
behaviour of the coefficients af^ one finds that the series has A < 00 if none of the solutions \ij of the 
characteristic equation (45) satisfies the condition 

< Imj/m^"^ - 1| < 1 , i = l,...,i?. (102) 
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If A = oo the scries is everywhere divergent and the expansion in factorial series is only formal; in this 
case another method must be used to calculate U'^"\l) (see section 8 on the Laplace's transformation). 

If A < oo the convergence of the factorial series is logarithmic, that is, if Sm{x) is the sum of the 
first m, terms of Eq.(lOl), one finds that |5„i(.x) — ^ ■m^~^ for large m,. The scries converges 

if a; > A, more and more quickly as x increases. As usually one finds A ~ 1, it is not convenient or 
possible to calculate J7(")(l) directly by summing the factorial series. Therefore, chosen a number x^ax 
conveniently large, one calculates C/'^"^(.t) for i? contiguous values of.T, U'^°'\x„iax), U''°'\x,y-,ax + '^), ■■■ , 
C/^"Ha;maa;+fi-l), where the series converges faster, and one uses repeatedly the corresponding recurrence 
relation, Eq.(29) or Eq.(31), in order to obtain the values of U'^°'\x) for x = Xmax — ^iXmax — 2, . . . up 
to a; = 1. A drawback of this procedure is that the recurrence relation may be unstable, so that each 
iteration causes a loss of precision. 

Let = minj The recurrence relation is unstable if 

1. A> 1: in this case each iteration increases the error on U'^'^\x) of a factor A. 

2. A=l, and /z^") is a root of Eq.(45) of multiplicity m > 1: in this case n iterations of the recurrence 
relation increase the error on J7(")(a;) of a factor n™-^; this is a kind of instability weaker than the 
preceding one. 

If the recurrence relation is stable Xmax can be chosen large at will, with no effect on the precision 
of [/(")(!). In the case of instability with A > 1, in order to obtain the result [/("^(l) with a number E 
of exact digits, the calculations of U'^°'\xmax + i) must be performed with a greater number of decimal 
digits C = E + Xjnax logio A. Supposing ai"^ ~ s!, a rough estimate of the number s^ax of terms of the 
series needed to obtain the sum with such a precision is Sfjiax A '^'-^ a^max* Performing the calculations 
with fixed precision arithmetic with C digits, it is convenient to choose a value of Xmax as low as possible 
in order to obtain the greatest E, compatibly with the rapid increases of Smax and of the computing 
time (see an example in section 7.2.4). Performing calculations with multiprecision arithmetic, C can be 
chosen at will, and Xmax can be increased; a convenient choice which minimizes the estimate of Smax is 
C/E 1 +lnA (see an example in section 9.4). Fixed C/E, by varying E one sees that Xmax oc E and 
Smax oc E. Therefore the number of terms of the series (and even the computation time) is proportional 
to the number of digits of precision of the results; this is true even in the stable case provided that one 

chooses Xmax OC E. 

6.2 Truncated expansion in e 

We are interested in the results in the limit D ^ A; therefore, defined e = (4 — D)/2, we expand all 
the quantities in e, truncating the expansions at the first terms, and we perform all the calculations 
using truncated series; in this way the coefficients of all the powers of e are found numerically, including 
the negative powers. This technique is perhaps not the most efficient, but is very versatile. We have 
implemented in the program SYS (see section 11) the arithmetic of truncated series, so that the use of 
series becomes as simple as with ordinary numbers. 

The time of computation of a factorial series grows approximately as n^, and it is mainly due to 
multiplications. The division between series is an operation less frequent than the multiplication as 
it occurs only once in the calculation of each a,, with the recurrence relations (48) or (57) or in the 
calculation of each U{x) using the recurrence relations (29) or (31). As effect of cancellations of terms 
when series are summed, and of divisions by series beginning with a non-zero power of e, first or last 
terms of the expansions in e may be lost; in this case the calculations umst be necessarily performed 
with an initial number rie of terms of the expansions greater than the desired final number n'^ of terms of 
U{x). The number is chosen empirically. The loss of terms of the expansions in e is frequent when the 
recurrence relations (29) or (31) are used to obtain U{x) for x < A; this is due to the fact that for such 
values of x the expansion in e of po{x) begins with a non-zero power of e. Furthermore, we found that 
it is critical to recognize the cancellation of the first coefficient of a series used as divisor when, because 
of the numerical errors, the coefficient is a very small number instead of zero; in this case a numerical 
cutoff must be carefully used. 
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7 Applications to simple one-loop integrals 

In this section we discuss in detail (as we assume the reader unfamiliar with these operator techniques) 
the solution of the difference equations for the one-loop vacuum and self-energy integrals. 

7.1 One- loop vacuum integral 

Defining 

^ ' J [k^ + mlY ^ ' 

we want to calculate the master integral J(l). According the single identity of the set [l; q] J{x) satisfies 
the homogeneous difference equation 

ml{x - l)J{x) -{x-l- D/2)J{x - 1) = . (104) 

We look for a solution^^ of this equation in the form of a factorial series (see section 4.3), 

oo 



s=0 



Introducing the operators tt and p, multiplying by x and using Eq.(39) we find the first canonical form 
of the difference equation 

{{uml - l)p2 + ((2/im? - l)(7r + D/2)p + iimlTz{-K - 1)) V{x) = . (106) 

The characteristic equation Eq.(45), /2(m) = 0, gives the value = l/m\. Substituting this value of ^ in 
Eq.(106) one obtains 

((tt - 1 + i:»/2)p + 7r(7r - 1)) V{x) = . (107) 

The indicial equation Eq.(49), /i(7r = 1 -(- JsT) = 0, gives the value K = —D/2. Using the recurrence 
relation (48) one finds 

_ T{s + D/2)T{s + Dl2 + l) 

- r{D/2)T{D/2 + i)r{s + if° ■ ^^^^> 

Prom Eq.(71) one finds that ao = (mf )^/^. The final result is 

J(x) - f - r(. + D/2n.s + D/2 + 1) T{x + 1) 

- y^^i) r(£>/2)r(£>/2 + i)r(s + 1) r(a; + 1 + D/2 + s) • ^ ' 

The coefficients behave for large s as Og/s! oc s^~^, therefore the term of the factorial series behaves 
as s^/^~^~^, so that the scries has abscissa of convergence A — D/2. This value signals the divergence of 
the integral (103). The recurrence relation (104), being of order one, is numerically stable. Therefore it is 
possible to obtain J{x) for a large integer x ^ D/2hy summing a few terms of the factorial series and to 
use the recurrence relation to obtain J(l). For a numerical example in the limit D ^ 4 see section 7.2.4. 

7.2 One-loop self-energy integral 

Let us consider the one-loop self-energy integral 

where Di = k'^ + ml and D2 = {p— kY + m\. A simple inspection of the system of identities [l . . . 2; 
shows that there are three master integrals, J[d^k\/ D1D2, J[d^k]/D\ and J[d^k]/D2. We must find 
three difference equations for 

^^The exact solution J{x) = aoT(x — D/2)/V{x) does not have the form of Eq.(33). 
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The solution of the system made up of the set of identities [l . . . 2; q° ^] gives the system of difference 

equations 

{x - D)I(x - 2) + {-p^ + ml - m\){2x -D- l)I{x - 1) + R^{p^, -ml, -ml){x - l)I{x) 

+ J(x - l){{D/2 - l){p^ +ml+ ml) - (x - 2){p^ + ml))/ml = , (112) 

ml{x - l)J{x) - (x - 1 - D/2)J{x - 1) = , (113) 

ml{x - l)K{x) - (x - 1 - D/2)K{x - 1) = , (114) 

where R was defined in Eq.(79). The difference equation for J(x) and K{x) has been solved in the 

preceding section, so that we consider here only the second-order difference equation for I{x). We discuss 
separately the massive case and the cases where one of the masses m\ or is zero. 

7.2.1 Case 1: rm 7^ 0, m2 ^ 

Following section 4.3 we substitute /(x) = ^j.^V{x) and J(x) = ^:''W(x) in Eq.(112), multiply it by 
x(x — 1) and insert the operators tt and p. We obtain the canonical form of the equation 

(/3(m)p' + /2(7r, ^x)p' + /i(7r, ii)p + /o(7r, V{x) = (53 (Tr)^^ + 52 (Tr)^^ + 5i(7r)p) W{x) ; 

(115) 

fi and Qi are polynomials in tt, h, D, , ml and m|, whose explicit expressions are not shown for brevity. 
The characteristic equation /s (/x) = has the two different roots 

M = M± = , . . \2 I 2 (116) 

(p ± im2Y +m{ 

(we define p — \pp^ , and if < then p — iyj—p^). The solution of the difference equation (112) is 

J(x) = 77+/^°(^) + ??-/5°(x) + /^^(x) , (117) 

where are the solutions of the homogeneous equation corresponding to /i± and is one solution 
of the nonhomogeneous equation. 

Let us now consider the homogeneous solution I^^ . We write I^'^{x) = n'^V^^{x) and we look for 
the coefficients of the expansion in factorial series of V_^'^ (x) 

00 

K_^°(a;)=^a,-p^-V (118) 

Substituting /x = the homogeneous part of the canonical form (115) becomes 

(/2"(^)P' + /r(^)P + /o"(^)) Vl'^^ix) = , (119) 

/2-(7r) =2im2p(27r + D - 5) , 

/r(7r) = ((p' + - mDiir + D - 3) + 2im2p(37r - 4)) (tt - 1) , (120) 
/o~(^) =(P^ + mi - 7712 + 2im2p)TT{TT - if ; 

the indicial equation f2 {K- + 2) = gives K- = (1 — D)/2. Fixed Oq = 1, the other coefficients a~ 
can be found using the recurrence relation (48); the behaviour of for large s can be determined by 
considering the recurrence relation as a difference equation in s for a~ and solving it; one finds 

a:/s\ « + C2B'/s , (121) 

where B = /x_/(/x_ — while Ci and C2 are constants. If \B\ > 1 the series (118) never converges 
(in fact the condition (102) is satisfied); if \B\ < 1, the series converges with abscissa of convergence 
X = D — 2. The solution I^'-^ can be obtained in analogous way with the replacement im2P —im2p. 
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Let us now consider tlie nonliomogeneous solution I^^{x). We write I^^{x) = fi'^V^^ (x). The 
value of 12 must be taken from the expansion in factorial series of J{x) obtained in section 7.1: /x = 1/mf. 
Substituting this value of n in the canonical form, Eq.(115) becomes 

(/3(^)p' + /2(7r)p2 + /i(7r)p + /o(7r)) F^^(x) = {m{7r)p^ + g2{^)p^ + 9ii^)p) W{x) , 

^ ' (122) 



/3(7r 

ff2(7r 
where J{x 



= (3(p2 + mlf + 2m\{p^ - m^)) tt - 5(p2 



+ (D-5)m2(/_^2) 



i?(p' - ml)) (Stt - 4) + TO2(p2 ^ ^2 _ ^2)(^ + £, - 3)] (tt - 1) , 



=i?2(p^-m?,-mi)7r(7r-l)2 , 



= + m^)(27r - 3) - {D/2){p^ +m\+ ml) 



(123) 



= (p2 + ^2)(^ _ _ (p/2)(p2 + ^2 ^ ^2) _ ^2^ (tT - 1) , 



(rni)~'^VF(x), W{x) = '}Z^=o^sP^^~'^ , Kb = —D/2, and the coefficients hg are the coeffi- 
cients as defined in Eq.(108). 

The right-hand side of Eq.(122) can be written as X^^g^s/^" 
expression of Cg is given by 



n.Ka-S 



, where Kc = K^, + 3, and the 



Cs = bsQsiKc - s) + bs-ig2{Kc - s) + bs-2gi{Kc - s) . 



(124) 



Then, the coefficients a., of the expansion of V^^{x) = CLsP^^ * can be found by using the 

recurrence relation (57); the large-s behaviour is 



(125) 



NH 



where B± = mj~^/(m^^ — /i±) and d are constants. If > 1 or \B-\ > 1 the expansion of V 
never converges (the condition (102) is satisfied); if \B+\ < 1 and < 1, the series converges, and the 
abscissa of convergence is A = max(£' — 2, D/2 — 1) if Ci and C2 ^ 0. 

Now we must determine the constants ry+ and ry_ ; we compare the large-a; behaviours of the solutions 
If o and I^", 



tHO 

-'± 

jNH 



{X) 

(x) 



(p2 + mi)-i(m2)^/2-x-^/2 p2 ^ _^2 ^ 
(2 - £))(4mi)-i(m?)^/2-^a;-^/2 = _^2 ^ 



(126) 

(127) 



obtained from the first term of factorial series, with the large- a; behaviour of I{x). We point out that 
these behaviours and the values of r]± later inferred are rigorously valid only if all the factorial series have 
finite abscissa of convergence; however the deduction may be extended to the cases where the series have 
A = 00 with the help of the integral representations of section 8.6.2. 

If p2 > —ml we are below the deformation threshold, and we can write I{x) as radial integral (see 
Eq.(80)) 



I{x) 



1 



rffc2(fc2)IV2-l^ 

r(D/2) J (k^+ml)^ pk 



the large-a; behaviour is given by Eq.(71) 
I{x) « (mf)^/2-^a;-^/2(p2 + m2)-i . 



(128) 



(129) 



The comparison of Eq.(129) with Eqs.(126)-(127) shows that, as fi± 7^ l/ni\, I^^ do not contribute to 
/, and that only I^^ contributes, so that 



??+ = = {p^ > -ml, mi ^0, 7712 7^ 0) . 



(130) 
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At the deformation threshold = —m^ one finds from Eq.(93) for large x 

I{x) « (m?)^/2-V2-a=a;(i-i?)/2^ . (131) 

2m2 

as /i_ = 1/mf and = —D/2 + 1/2, Eq.(98) is satisfied, so that I^*-* contributes to /, and rj- is 
different from zero. Comparing Eq.(131) and Eqs.(126)-(127) one finds 

r?_ = mf "^m^^A/7r/2 , /?+ = [p"^ = -ml, mii^Q, m2T^Q) ■ (132) 

If < —m^ we are above the deformation threshold, and the radial path of integration turns around 
the point k"^ = {p — im2Y- We split I{x) = Ia{x) + h{x) (see Eq.(81)), where is given by Eq.(128) 
and lb is 

1 / (133) 



r(D/2) j {k^ + mff \pk ^ ' pk 

(p— im2)^ 

Noting that for large x the contribution to this integral comes from the neighbourhood of the singu- 
larity, one finds 

h{x) « ^(i-^)/2^ ( !!^') , (134) 

m2 V / 

so that 

r?- = ^(^^ , V+ = (/<-mi mi^O, m2^0). (135) 

7.2.2 Case 2: mi 7^ 0, m2 = 

In this case we have 11+ = H- = + m^). The second-order difference equation (112) simplifies to 

the first-order equation 

{D-x~ l)I{x - 1) + (x - l)(p2 + ml)I{x) -{x- l)J{x) = , (136) 
whose solution is of the kind 

I{x)='nl"°{x)+I''"{x) . (137) 
In a manner analogous to the previous case, one obtains the expansion in factorial series 

oc 

7«0(^) ^ (p2 ^ ^2)-. ^ a^p2-D-. . (13g) 
s=0 

the large-s behaviour of the coefficients is a^/s! oc s^^~^, so that the series has abscissa of convergence 
X = D — 2. Considering the nonhomogeneous solution, 

00 

/^^(x) = (m?)-^a,p-^/2-% (139) 

s=0 

the coefficients have the largc-s behaviour 

a,/s! ^ Cis(3^-6)/2 + C2S^-2 + , (140) 

where _Bo = 1 + m\lp^ and are constants. If > —m\j2, \Bq\ > 1 and the series does not converge; 
if P^ < —ml/ 2 the series converges with abscissa of convergence A = max(£) — 2,£'/2— l),ifCi7^0 and 
C2 7^0. 
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Now wc determine the constant rj. The large-a; behaviour of the solutions is given by the first term 

of the factorial scries 



I"'^{x) « (p2 + m2)-^2;2-^ , (141) 

{„-2( 2\D/2-x -D/2,f „2 n 
p (m^) X +. p > U , 

In the whole region > the factorial scries expansion of never converges so that the large-a; 
behaviour shown above is valid only in asymptotic sense. Using the integral representations of sec- 
tion 8.6.2 one can show that if > the large-a; behaviour of contains an additional contribution, 
exponentially small in comparison with the main contribution 

7^^ (a;) « p-2(m?)^/2-a;-^/2 + T{D/2 - 1)( + m2)^-2-a;2-^. (143) 

The large-a; behaviour of I{x) can be deduced from Eq.(128) and Eq.(133), as in the massive case. 
But if > there is a difference: the integrand of Eq.(128) is discontinuous in the point fc^ = p^, in fact 

p^ + fc^- ^ \k/p k<p, 

2pk \p/k k>p. ^ ' 

The presence of this discontinuity gives rise to an exponentially small additional contribution proportional 
to [p^ + m\)~^ in the large-a; behaviour of I{x). So we split the integral over k"^ into two parts, 

~ T{D/2) J (fc2 + m?)^ p2 ^ 



^T{D/2)J (P+m?)^ Vfc^ V^V rWJJ' 

each one having a different large-x behaviour. One finds that the large-a; behaviour of I{x) is identical 
to that of I^"{x) for p^ > 0, Eq.(143), so that 

,7 = (p2 > 0, mi ^ 0, m2 = 0) . (146) 

If p^ < the large- a; behaviour of I{x) can obtained from Eq.(133) (with m2 = 0). One finds the 
same result found for p^ > 0. Therefore for p^ < the exponentially small contribution is present in the 
large-a; behaviour of I{x) but not in that of I^^{x): it must come from 'ql^'~'{x). Therefore the constant 
?7 is 

r, = V{D/2 - 1) ((p2 + mf)/(-i p)f~^ (p^ < 0, rm 7^ 0, ma = 0) . (147) 
7.2.3 Case 3: mi = 0, m2 

In this case the denominator raised to x has zero mass, and the term containing J in Eq.(112) disappears, 
so that the difference equation becomes homogeneous; the solution is of the kind 

/(.x) = r,+/:^0(x) + . (148) 

I^'^ix) are the same functions considered in section 7.2.1 with rrii = 0, and their large-x behaviour was 
given in Eq.(126). But the large- .t behaviour of I{x) cannot bo found as in the massive case; in fact if 
mi = the integrals (128) and (133) arc strongly divergent in fc^ = for large x and only dimensional 
rcgularization can give a finite value to I{x). The large-a; behaviour of I{x) must be determined by other 
methods, for example by writing I{x) as integral over one Feynman parameter and using the saddle-point 
method; one finds 

= o ^ n/o^ (T^PM±/m2)(^-^^/' (mi = 0, m^ ^ 0) . (149) 

2m2 sm(7rL'/2) 



28 



X 


J(x) 


I^^(x) 






9 


017857 + 033442f 


047713 -\- 089160f 


—0 008928 


— 007323f 


8 


023809 + 040621^ 


059006 + 100337f 


—0 011904 


— 007663f 


7 


n 033333 + 050203f 


075609 + 11 3249f 


—0 016666 


— 0071 59f 


6 


05 + 062805f 


101857 + 126598? 


—0.025 


— 003929f 


5 


083333 + 076898e 


148085 + 133074e 


-0 041666 


+ 009010e 


4 


0.166666 + 0.070464e 


0.245635 + 0.090010e 


-0.083333 


+ 0.067207e 


3 


0.5 - 0.288607e 


0.548843 - 0.442122e 


-0.25 


+ 0.534955e 


2 


- 0.577215 


0.282094e-i + 0.519388 


-0.25e-i 


+ 0.144303 


1 


-e-i - 0.422784 


0.282094e-i - 1.645293 


0.75e-i 


+ 1.067088 








-0.564189e-i - 0.238530 


-0.5e-i 


- 0.211392 



Table 1: Values of J{x), I^°{x) and I^"{x). 



7.2.4 Numerical example 

Let us describe in some detail the calculation of /(I), J(l) and K{1) in the case —p^ = m\ = rri^ = 1- 

As mi = 1712 then K{x) = J{x), and wc consider only J{x) and I{x). The root of characteristic equation 
of the difference equation for J{x) (Eq.(113)) is /z = 1; the roots of the characteristic equation of the 
homogeneous part of the equation for I{x) (Eq.(112)) are ^_ = 1 and /U+ = —1/3. I{x) is written as 
sum of the homogeneous and nonhomogeneous part, Eq.(117); from Eq.(132) one finds 77+ = 0, so that 
I^'^ does not contribute. Therefore we have to evaluate three functions: J{x), I^'~'(x) and I^^{x) for 
X = 1 and £> — > 4. The factorial series expansions of all the three functions have abscissa of convergence 
A ^ 2 for — !■ 4 and do not converge for x = 1. Following the procedure of section 6 we evaluate the 
series for a large value x = Xmax and, as the equation for I{x) is of second order, even for x = Xmax + 1- 
The convergence becomes faster by increasing Xmax, but, a.s A = = 3 > 1, the recurrence 

relation for I{x) Eq.(112), rewritten in order to obtain I{x — 2) from I{x — 1) and /(x), is unstable. 
Each application of the recurrence relation increases the error of the value of I{x) of about a factor 3. 
Therefore we must choose a value of Xmax which is a compromise between speed of computation and loss 
of precision in the result. We choose Xmax = 8. Calculations are performed by setting D = 4 — 2e and 
expanding in e, truncating the expansions at the first three terms. Doing the calculations with 19 digits, 
the convergence is attained for all the factorial series for a; = 8 (and consequently for a; = 9) in about 
4000 terms. The values for x < 8 are calculated by using repeatedly the recurrence relations (112) and 
(113). The application of the unstable recurrence relation enlarges the error of about a factor 3^ « 2000, 
corresponding to a loss of about 3 significant digits in /(I). Values of J(x), I^^{x) and I^^{x) are 
shown in Table 1, with only the first two terms of the expansion in e, and the coefficients with only 
6 digits to save space. The value of J(l) calculated agrees with the exact result — + 7 — 1 + 0{e) 
within the precision of the calculation (7 is the Eulcr's constant). Inserting and the 

value ?7_ = ^/n/2 (from Eq.(132)) in Eq.(117) one finds 

7(1) = (1 - 2 X 10-^^)e-^ - 0.3910150291357503+ 0(e) , (150) 

with an error of 4xl0~^^ on the constant term in comparison with the exact result e~^+2— 7r/\/3— 7+0(e). 
A numerical value of r/_ can be obtained independently, using the identity (100) 

7(0) = + 7^-^(0) ; (151) 

inserting the values of 7(0) = 7s:(l) = J(l), l"O{0) and 7^-f^(0) one finds 

r]- = 0.8862269254527578 + 5 x lO'^^e + 0{e^) , (152) 

with an error of 2 x lO""*^^ on the constant term. Using this numerical value in the place of the analytical 
value one obtains a value of 7(1) with the same precision as Eq.(150). In Table 2 wc show for different 
choices of x^ax the number of terms of series necessary to evaluate J(8), 75*^(8) and I^"{8) with 19 
digits of precision and the obtained values of the finite part of 7(1); it is evident that by increasing x^ax 
the series converges faster, but the precision degrades because of the increasing number of applications 
of the unstable recurrence relation. 
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^raax tCrinS 



finite part of /(I) 



30 
25 
20 
15 
10 
9 
8 
7 
6 



125 
154 
217 
395 
1470 
2454 
4439 
13086 
36210 



0.3910008887063124 
0.3910149952724784 
0.3910150292106927 
0.3910150291388126 
0.3910150291357554 
0.3910150291357472 
0.3910150291357503 
0.3910150291357507 
0.3910150291357507 



Table 2: Dependence of the finite part of /(I) on x. 



8 Solutions of difference equations by means of Laplace's trans- 
formation 

The expansion in factorial series is certainly the most direct method of solution of difference equations 
with polynomial coefBcients; however, for some values of masses and external momenta of the diagram, as 
we have seen in the above example, the abscissa of convergence of factorial series may become infinite. In 
this case the factorial series become divergent for every value of x and therefore useless, so that another 
method of solution must be used: the Laplace's transformation method[19]. 

This method is described in section 8.1, and applied to the systems of difference equations in section 
8.2, 8.3 and 8.4. Techniques used for integrating the differential equations obtained from the application 
of the method are described in section 8.5. The application to simple one-loop integrals is shown in 
section 8.6. 

8.1 Transformation of a difference equation 

Let us consider the difference equation 



where Pi{x) are polynomials in x of maximum degree P. The Laplace's transformation method consists 
in the substitution 



where I is a line of integration suitably determined and where v{t) is found from a certain differential 
equation. Writing the coefficients as 



Pq{x)U{x) +pi{x)U{x + 1) + . . .+pn{x)U{x + N) = , 



(153) 




(154) 



p 




(155) 



substituting Eq.(154) in Eq.(153) and integrating by parts one finds 



N . P 



Y,Pk{xMx + k)= / dtt--i^$,(0(-t)V*)(t) + [7(a;,i)], 



(156) 




where 



N 




(157) 



fe=0 




j=0 m=0 



(158) 
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Eq.(154) provides a solution of the difference equation (153) if v{t) is a solution of the differential 
equation 



J2^i{t){-tyv^'Ht) = , (159) 

provided that the line of integration I be chosen so that I{x, t) has the same value at each endpoint of the 
line, if the line is open. Note that the difference equation (153) has order TV with coefficients of degree 
P, while the differential equation (159) has order P with coefficients of degree N + P. 

The singular points of this differential equation are 0, oo and the zeros ti (of multiplicity nii) of the 
characteristic equation 

^p{t) = ; (160) 

these points turn to be always regular singular points in the case of the differential equations encountered 
in this work. 

We choose as lines of integration the lines which begin in the origin and end in one of the singular 
points ti. This is a convenient choice. One can show that I{x, t) = at each endpoint of such lines, if the 
integral over t of Eq.(154) is finite. Under these conditions, U{x) is a solution of the difference equation. 

We can construct a set of rrii = N functions Uij{x) which form a fundamental system of solutions 
of the difference equation (153) by defining 

Uij{x) = j dt f'-^Vijit) 3 = 1,... , mi, (161) 

li 

where Vij{t) is one of the mi solutions of the differential equation singular in ti, and U is the line which 
begins mt = Q and ends mt = ti. 

It is important to note that the characteristic equation (160) of the differential equation and the 
characteristic equation of the difference equation (153) turn out to be identical, so that we can readily 
identify the singular points ti of the differential equation with the solutions yUj of Eq.(45). 

Now we consider the nonhomogeneous equation 

Af N' 

J2Pk{3:)U{x + k)=J2^k{x)T{x + k) , (162) 

k=0 fe=0 

where qkix) are polynomials and T{x) is a solution of some difference equation. A particular solution 
U^^{x) is found by substituting into the equation 

T{x) = j dtf'-^'wit) , U^"{x) = j dtt^'-'^VNHit) , (163) 

It It 

where It is a known line of integration and w{t) is a known function, solution of a differential equation 
analogous to Eq.(159), obtained from the difference equation satisfied by T{x). Provided that I{x, t) = 
at the endpoints of It, one finds the nonhomogeneous differential equation 

Y.MtK-tyv'i^^Hit) = j2'^im-tyw^^\t) , (m 

i=0 i=0 

whose solution gives Vnh- Then U^^ is found using Eq.(163). 



8.2 Transformation of the system of difference equations 

In section 3.2 we described the construction of the system of difference equations. The algorithm used 
yields a system in triangular form; applying to this system the Laplace's transformation one obtains a 
system of differential equations with the same triangular structure and the same ease of solution. There 
is, however, a complication: the Ith difference equation 

^Pk{x)Ui{x + k) = J2Yl l3k{x)Uj{x + k) (165) 

fc=0 j=0 fe=0 
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of order N and with polynomial coefficients of degree P is transformed, using the Laplace's transformation, 
into a differential equation of order P 



p i-i Pj 

J2Mt)i-tyvl'\t) = W (166) 

i=0 j=0 i=0 

and coefficients of degree N + P. As effect of the particular algorithm of construction and solution of the 
system of identities, N is usually small (typically 1^4) while P may be large (3 ^ 30). Therefore the 
differential equation may have a high order. Calculations in some test cases have shown that solution of 
high order equations slows down the calculations and is source of undesired numerical errors, so that, if 
possible, it is better to avoid it. 

We have discovered that this difficulty can be overcome by applying the Laplace's transformation to 
the identities obtained by integration-by-parts before the insertion in the system of identities, building a 
system of "transformed" identities between the "transformed" integrals v{t) instead of the real integrals 
U {x) . The solution of the system of transformed identities will provide a system of differential equations 
of smaller order; a small price to pay is the (possible) appearance of spurious singular points in the 
equations so obtained (see section 8.5.2). More in detail, a generic integral is transformed into 

f TT^p-" / . ^ irred f 

Uniafsix) = / [d^fcl] . . . [d'^kN,] = / dt t^-'Vni^0{t) , 

ii h ' ' ' in "J 

where the line I is unspecified; the values of the functions I{x, t) at the endpoints of the line are always 
zero because of dimensional regularization. Therefore, a generic integration-by-parts identity, 

niaP 

{rniai3 and Sniafi are independent of x) becomes the transformed identity 

({Sniafj - rnia0(Xl)Vnia/3{t) - t r„ic»/3 ^^^^"^ (^) j t"^ = (167) 

which is a differential equation between the functions Vniafiit)- 



8.3 Construction of the system of differential equations 

We want to build a triangular system of differential equations 

Pi i-i Pik 

T.Piiit>^\it) = Y.T.'lMt)vlil{t) , 1 = 1,. ..,L'^, (168) 

i=0 fc=l j=0 

between a set of "master transformed" functions Vmi(t) (to = 1, . . . , Nd — Nk + 1) analogous to the system 
of difference equations (25) between the master functions Umi{x) discussed in section 3.2. For this reason 
we devise an algorithm of construction and solution of the system of identities similar to the algorithm 2: 

Algorithm 3 Consider the algorithm 2 with the following modifications: 

1. After the step 8 of algorithm 1, transform each integration-by-parts identity using Eq.(167); in 
the subsequent steps replace everywhere the integrals /[rf^fci] . . . [d^kNk]Wniai} with the functions 

Vnia/Sit)- 

2. The transformed identity obtained does not depend on the index ai ofW in the step 7 of algorithm 1 
because it appears as an overall factor ; in order to restore the total number of different identities 
we have decided (somewhat arbitrarily) to differentiate ai times with respect to t each transformed 
identity. 

3. Ignore step 4 of algorithm 2. 
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4- Derivatives of master transformed functions have priority of extraction lower than other generic 

transformed integrals. 

5. Add the new entry the greatest derivative to the list of priorities after the entry 9(h)vi of algorithm 1. 

With a suitable choice of the parameters Oi and bi (see the end of section 3.2), by means of this 
algorithm we can identify the master transformed functions Vmi (t) as the functions which satisfy equations 
of non-zero order, and we can work out a set of differential equations among them. If each function 
Vmi{t) corresponds to an integral with a different combination of denominators, the system is obtained 
directly in the triangular form (168); if, on the contrary, there are different master transformed functions 
Vm^i+iit), ■ ■ ■ ,'Vm,i+G{t) Corresponding to integrals containing the same combination of denominators, the 
algorithm provides a set of G simultaneous differential equations containing all these G functions, which 
are conveniently transformed into triangular form using a procedure quite analogous to that described in 
section 3.2. As final result one obtains a system of differential equations with the triangular form (168). 
It is important to note that the functions 

Fml{x) = J dtt^-^Vml{t) (169) 

I 

arc not necessarily identical to the master functions U„ii{x) defined in Eq.(22). While v„ii{t) is a master 
transformed function and satisfies a differential equation of non-zero order, Fmi{x) satisfies a difference 
equation of different order, which may be zero; if so, Fmi{x) is not a master function. This fact frequently 
occurs when many master integrals have the same combination of denominators. For these reasons the 
number and the structure of the master transformed functions Vmi{t) must be found independently of 
the master functions Umi{x). See section 9.4 for an example. 



8.4 Correspondence with factorial series and initial conditions 

Let U{x) = Yl^Uaix) be the solution of a generic difference equation; the initial conditions for the 
integration of the differential equation can be determined by comparing the integral representation of 
Ua{x) with the factorial series expansion 

U^{x) = dt t^-W{t) = Ma E "-P"^""' • (170) 
The behaviour of Va it) near the singular point t = fia 

Va{t) W Aoail^a - t)^^ , t ^ fJLcc, (171) 

can be deduced from the known behaviour of Ua{x) for large x. Substituting Eq.(171) in Eq.(170), 
integrating over t and comparing the result with the large-a; leading behaviour of the first term of the 
series Ua{x) « ^'^aoaX^" one finds JsT^ and Aoa- 

= - 1 . ^Oa = aoa/{t^a"T{K^ + 1)) . (172) 

If necessary, the subsequent coefficients of the series expansion of Va{t) can be deduced in the same way 
from the coefficients asa- 



8.5 Integrating differential equations 

The choice of a effective numerical method of integration of the differential equations obtained by applying 
Laplace's transformation is not simple. In fact we must consider that: 

• The initial and final point of the path of integration / arc singular points of the solutions; the 
numerical methods usually used to integrate differential equations (for example, the Runge-Kutta 
method) cannot be used in singular points. 

• Master functions, coefficients of differential equations and (sometimes) singular points depend on 
D and are represented by truncated series in e. 
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• The method must be able to provide very high-precision vahies, and the time of computation must 
grow Unearly with the number of exact digits of the result; by using fixed-order methods (hke the 
Runge-Kutta method) the time grows exponentially. 

Therefore we have decided to solve the differential equations by using power series, expanding the general 
solution around a number of selected points near or on the path I and equating the sums of the scries in 
some intermediate points. Each power series will be evaluated inside the respective circle of convergence, 
so that the number of terms of the series necessary to attain a precision of E digits in the results (and 
also the computation time) will be proportional to E. 

8.5.1 Integration over the path 

The line of integration can have any shape, but for convenience it is assumed here to be the segment 
[0,/u]. The segment is subdivided into M intervals [tM+i,tM], [tM ,tM-i], ■ ■ ■ , [^2,^1], where ti = fj, and 
tM+i = 0; the (possible) singular points placed on the segment are = tp < . . . < tf < tf < /j. The 
choice of a line of integration which passes through singular points of the differential equation, instead 
of avoiding them, may be convenient: it allows one to check if v{t) is regular or singular in these points 
and it speeds up the calculations, avoiding the use of complex numbers. Let us consider the first interval 
[^2,^1]- The solution is expanded around the point t\ = n (which may be a singular point). The first 
coefficients of the expansions of the P solutions corresponding to the roots of the indicial equation arc 
obtained using Eq.(172). All the subsequent coefScients are obtained using recurrence relations obtained 
by substituting the expansions into the differential equation. Then the power series are evaluated in the 
point t2 = ti — ri/2, where ri is the radius of convergence of the scries and the factor 1/2 has been chosen 
in order to minimize the total time of calculations. In the next interval the solution is expanded around 
the regular point t2, the first P coefficients of the expansion are obtained from the already known values 
of V and its derivatives in t2, the subsequent coefficients arc obtained using the recurrence relations, and 
the power series is evaluated in a point = t^ — r2/2. This process involving expansions around regular 
points continues for m steps until a point tm is reached, placed inside the circle of convergence of the series 
with center the singular point tf, at a distance less than one half of the radius rf of the circle. Now we 
expand the solution around the singular point tf. We write the general solution as v{t) = CjVj{t), 
where each Vj corresponds to one of the P solutions of the indicial equation. The values of Cj are found 

by solving the linear system^"^ 'l2f=i '^j^j^\tm) = v'^^Him), i = 1,... ,P- The right-hand sides contain 
the already known values of v and its derivatives in 1,^: the necessary values of partial solutions Vjit^) 
are worked out by summing the corresponding expansions about tf', whose coefficients are found using 
recurrence relations. Then the singular point is got over by evaluating v (and its derivatives) in the 
regular point tm+i = if — rfY2. The process is repeated until the next singular point tf is reached, then 
tf, etc., up to the final point tp = tM+i = 0. 

The integration of v{t) over t needed to obtain U{x) can be easily carried out by integrating the 
expansions in series in the corresponding intervals: 

^ M *'+! P 00 

U{x) = dt t'^-^vit) =^ dt t^-i ^^a(*'^)(t - ti)^''+' . (173) 
The integrals 

i+b b 

s)= J dt t^{t - t)>'+' = Jdy{y + t)iy''+' (174) 

i+a a 

which appear in Eq.(173) can be expressed in terms of incomplete Beta function (note that k is not an 
integer). If j = the integral is immediate; ii j is positive integer the value can be efficiently computed 
using the recurrence relation 

s) = I{j - 1, ,s + 1) + tl{j - 1, s) . (175) 



^ ^ It is important to note that if m + 1 solutions of the indicial equation coincide for D ^ 4, the Wronskian determinant 
of the linear system is proportional to {D — 4)'"; this causes a loss of m terms in the expansions in D — 4 of cj and, 
consequently, in the expansions of v in the rest of integration, and therefore in U{x). This mishap sometimes occurred in 
the final singulcir point t = in the calculations of section 9.2. 
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8.5.2 Singular points depending on D 

The coefficients ai^'"'^ and the exponents Kij in Eq.(173) depend on D; therefore all quantities are 
expanded around = 4, and truncated series are used in the calculation, as described in section 6.2. 
A new feature, characteristic of the differential equations obtained by solving the system of transformed 
identities, is the appearance of spurious apparent singular points t, not corresponding to any solution 
H of the characteristic equations of the difference equations obtained by solving the system of 'original' 
identities. These apparent singular points are depending on D, and correspond to regular points of the 
solution of differential equation, in contrast with regular singular points, which are independent of D. In 
general the line of integration can be deformed in order to avoid these spurious singular points; but if for 
£) — > 4 one or more of these mobile points tends to one of the endpoints of the line, t = n ov t = 0, we 
cannot avoid it, and we encounter difficulty in working out the sohition near these coalescing points. Let 
us explain this fact, by considering a homogeneous differential equation with polynomial coefficients 

P Oi 
i=0 j=0 

and supposing for simplicity that the equation has only one apparent singular point to{D) such that 

to{D) = 0{D - 4) for D ^ 4 . (177) 
The coefficients of the expansion v{t) = J2^o «si*^^ can be found using the recurrence relation 

_ j:T=i^s-,.f,{K+s-j) 

- MkT7) ' ^'"'^ 

where = if i < 0, m, = max^ gi, fj{k) — X^iLo ^(^ — 1) • • • (fc — i + , and K is one of the roots of 
indicial equation fo{K) = (note the analogy with the solution of a difference equation with expansions 
in factorial series). 

As the point to is a regular point, the coefficients a,,, which arc functions of D, have values in D = 4 
finite and, in general, different from zero; unfortunately, there are problems for calculating them. As a 
consequence of Eq.(177), the function /o vanishes if D = 4, and the other fj do not vanish (but the sum 
in the numerator of Eq.(178) always vanishes for D = 4). so that the recurrence relation (178) turns out 
to be very unstable, with a degree of instability proportional to 1/(1? — 4) . Performing the calculations 
of tta using truncated expansions in D — 4, each iteration of Eq.(178) (with increasing s) yields one new 
coefficient a^, whose expansion in £) — 4 has a number of terms reduced by one in comparison with as_i; 
after a few iterations, the given number of terms of the expansion in D — 4 is exhausted. We found that 
a solution of the problem is to modify the recurrence relation to 

~ h{K + s-l) ' 

so that the denominator is /i, which does not vanish for D = 4. The new recurrence relation requires as 
input the value of the sth coefficient before that its value is obtained; therefore Eq.(179) must be seen as 
part of an iterative process: 

1. set af^ = for s = 1, 2, . . . , Smax', 

2. apply Eq.(179) for s = 1,2,..., Smax obtaining the coefficients ai^^; 

3. repeat n times the step 2 until lai"^^-* — ai"''| = 0{{D — 4)™) for every s, where m is the desired 
number of term of the expansion in Z) — 4; the convergence is guaranteed in about m steps by the 
fact that /o = 0{D - 4). 

Analogous modifications must be made to the recurrence relation in the case of nonhomogeneous 
equations, or in case of two or more spurious singular points coalescing to the same endpoint. An 
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example of equation with one mobile singular point is the third-order equation^^ 

- 2{t - l){8t + 1)(4(£) -l)t-D + 4))t^v'^{^ + {-32{D - 1){5D - 7)t^ 

+ 8{D - 1)(25D - 63)^2 + {-4D^ + lOOD - 192)t - (D - 4)(9D - 22))t^V2^ 
+ (-32(D - 1){D - 2){3D - 5)t^ + 32(D - 1){D - 3)(5L' - 12)t^ 
+ {18D^ - 70D^ - AAD + 288)t + {D - A){-ViD^ + 12D - 100))to2h 

+ 2{D - A){D - 3)2(16(1? - l)f^ + 8Dt -30 + 10)?;2h = , (180) 

which is the homogeneous part of the equation satisfied by the function V2h{t), corresponding to the 
integral of Fig. 2h of section 9.2, used in the calculation of Eq.(208) (Eq.(180) remains the same for both 
choices of the line Di). The singular points are t = 1, —1/8, and {D — 4)/{4D — 4). The last singular 
point, which satisfies the condition (177), is a regular point with exponents 0, 1 and 3. The characteristic 
equation of the homogeneous part of the corresponding difference equation, 

{x-D + l){x -2D + A){2x - 3LI + 6)(3a; - AD + lQ)U2\,{x - 1) 
+ 2(a; - £> + 2) (21a;^ + (136 - QlD)x^ + (243 + QAD'^ - 253£))a; 

- 8(£> - 1){D - 2){2D - 5))i72h(a;) 

- 8x{x -D + 2,){2x - 3£> + 7)(3a; - AD + 7)C/2h(a; + 1) = , (181) 

has only the roots 1 and —1/8. 

8.6 Applications to simple one-loop integrals 

Now we consider the solution with Laplace's transformation of the difference equations analyzed in the 
examples of sections 7.1 and 7.2. 

8.6.1 One-loop vacuum integral 

Considering the integral J{x) of Eq.(103), the endpoints of the line of integration are the origin and the 
root of the characteristic equation fi = l/m\, so that we write 

J{x)= j dtf'-^jit) , (182) 



where vj satisfies the differential equation 

- t{mlt - l)vj{t) + {D/2 - mlt)vj{t) = . (183) 
The solution is 

vj{t) = C{l/ml - if'^-H-^'^ ; (184) 

the constant C can be deduced from the value of oq (see section 12), the value oi K = — 1?/2 and the 
relation (172). One finds 

C = (m?)^/2-i/p(£)/2) . (185) 

8.6.2 One-loop self-energy integral 

Here we consider I{x) of Eq.(lll), and in particular the case of non-zero masses. The homogeneous 
solution can be written using the Laplace's transformation as 

I^°{x) = J dte-^VHo{t) . (186) 





The equation (180) was found by solving the system of transformed identities. Applying the Laplace's transformation 
directly to the difference equation (181) one gets a higher order differential equation which docs not have the mobile singular 
point. This equation can be also derived from Eq.(180) by writing (3td'P{t)/dt + 4(D - 4)*(t))/(4(D - l)t - (D - 4)) = 
where *(t) is the left-hand side of Eq.(180). 
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The function VHo{t) satisfies the differential equation 

- t^i{t)v'Ho{t) + Mt)vHo{t) = , (187) 

where 

$i(i) = R\p\ -ml -ml){t - n+){t - M-) , 

^o{t) = -t'^R^{p'^,-ml,-ml) + {D-l){p'^ + ml-ml)t + 2-D . 
The solution of this equation is 

VHoit) = C ((m+ - t){^i. - t^-D _ (^gg^ 

Values of C such that are the same functions considered in section 7.2.1 are 

C± = 11^-^^'^ (mt - l^±f-''^'^ /T{{D - l)/2) . (190) 
Considering now the nonhomogeneous solution 

7^«(a;)= j dte-\NH{t) , (191) 



the function vnh satisfies the differential equation 

- t^i {t)v'j^H [t] + $0 {t)vNH [t] = t(t>x {t)v'j (t) - 00 {t)vj (t) , (192) 
where $o and $i are given in Eq.(188), and are 

Mt) = t{D{p^ +ml + ml)/{2ml) - 1) , 
and vj is given in Eq.(184). 



8.6.3 Numerical example 

Wc consider the calculation using Laplace's transformation of J(l), and numerically 

calculated in section 7.2.4. The corresponding differential equations arc Eq.(183), Eq.(187) and Eq.(192). 
The singular points of the system arc t = 1, —1/3 and 0. The equations arc integrated using the method 
described in section 8.5. The integral over t is divided into 4 intervals, with cndpoints 0, 1/8, 1/4, 1/2, 
1; a cutoff A <C 1 is conveniently introduced in the first interval [A, 1/8] because the solutions arc not 
regular in t = 0. Doing the calculations with 19 digits, the convergence of the expansions in f, is attained 
in about 80 terms; the finite values of J(x), I^'^(x) and I^^{x) for a; = 3 and :r = 4 are obtained by 
calculating the integrals with Eq.(173); the unstable recurrence relations are used only to calculate the 
values for x < 2, where the integrals are divergent, so that the error on /(I) is reduced by two orders of 
magnitude with respect to Eq.(150). 



9 Application to multi-loop diagrams 

After the self-energy diagram discussed in section 7.2, now we consider more complicated diagrams: the 

vacuum and the self-energy diagrams up to three loop, shown in Figs. 1 and 2, and the vertex and box 
diagrams up to two loops, shown in Figs. 3 and 4 (we have considered all diagrams such that the scalar 
integral (18) is always a master integral which does not factorize in a product of simpler master integrals). 
A complete discussion of all these diagrams would be rather long and we postpone it to future papers. 
However, to give an idea of the kind and complexities of the equations involved in the calculations, we 
will show some results. 

First of all. we distinct in each diagram g the topologically different lines which arc indicated in the 
figures with a number; of course in diagrams without numbers all lines are topologically equivalent. For 
each topologically different line I we set Di equal to the denominator of such line and we consider the 
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(a) (b) (c) (d) (e) 



Figure 1: Vacuum diagrams up to three loops. 



difference equation satisfied by tfie scalar master function Ugi{x) and the differential equation satisfied 
by the Laplace- transformed function Vgi{t), both corresponding to the scalar master integral, 



9.1 Arbitrary case 

In the left half of Table 3, considering arbitrary (non exceptional) values of masses and momenta, we 



• The diagram considered and, if present, the indication of the possible values of the index of the 
topologically different lines. 

• The number Uh of master integrals containing all the Nd denominators, determined with the pro- 
cedure of section 2.5, subdivided according to the number of scalar products in the numerator; for 
example, 1,4,3 means that we have found 8 master integrals altogether, of which 1 with no scalar 
product, 4 with 1 scalar product and 3 with a product of 2 scalar products. 

• In the Rc column we list the values of the order R of the difference equations in x satisfied by 
Ugi{x), for each possible choice of Di as one of the topologically distinct lines; the index C, where 
present, indicates (assuming values of external momenta below the deformation threshold) the 
number of constants rjj which are different from zero because the corresponding partial solutions of 
the homogeneous equation satisfy the condition (73). 

• The order S of the differential equations in t satisfied by the function Vgi{t), for each possible choice 



The orders shown within parenthesis are estimated from the subsystems of simultaneous equations 
(26), avoiding the transformation into triangular form; the index C is not shown in these cases. Analyzing 
the data of the table we observe that 

1. The differential equation has order less than the order of the difference equation; this is expected 
as Vgi{t) is an object simpler than Ugi{x). 

2. The order of the differential equation in t is equal to the number of master integrals, with the 
curious exception of the diagrams 2f3 and 2gl where S < nb- 

3. The choice of the line heavily affects the order and the complexity of the equations, as in the case 
of the diagram 2f, where the difference equation may have order 5,15 or 2. 

Due to the present limitations of the program used for the calculations, the values listed in this part of the table were 
calculated by giving arbitrary rational values to square masses and scalar products of external momenta. It is possible, but 
very unlikely, that the chosen values correspond to some particular case so that the results obtained do not correspond to 
the real arbitrary case. 




(194) 



lisfi^ 



of Di. 
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Figure 3: Vertex diagrams up to two loops. 




Figure 4: Box diagrams up to two loops. 
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Table 3: Number of master integrals and orders of the equations for the diagrams. 
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4. Some similarities appear between vertex diagrams, self-energy diagrams and vacuum diagrams with 

different number of loops 

(a) between a vertex diagram and the self-energy diagrams obtained by connecting two external 
vertices of the vertex diagram with a line and inserting an external momentum in this new 
line; 

(b) between a self-energy diagram and the vacuum diagram obtained by connecting its two external 
lines. 

The total number of master integrals and the number of master integrals subdivided according to 
the number of scalar products turn out to be the same, as well as the orders of the equations for 
the lines present in both diagrams. For example, compare the diagrams (3a, 2d, Ic), (2a, lb), 
(2b, Ic), (2c, Id), (3b, 2j, 21), (3c, 2o, 2p), (3d, 2r, 2s), (3e, 2t, 2u) and (3f, 2v). The similarities 
presumably exist even between three-loop self-energy diagrams and four-loop vacuum diagrams 
(some preliminary results seem to confirm it). Perhaps there is a relation with the heuristic "rule 
of the mapping" described in [1] for massless self-energy diagrams. 

5. The number of master integrals grows probably exponentially^^ with the number of loops and 
external vertices, and may be large. Consider for example the class of i-loop "sunset" self-energy 
diagrams with L + 1 denominators, shown for one. two and three loops respectively in Fig. 2a, 2b 
and 2e. The number ni){L) of master integrals for L = 1 to L = 5 is 1,4,11,26,57, respectively 
(the last two values come from preliminary analysis). These values seem to follow the law nb{L) = 
2^~^^ — L — 2 = J2i=2 i^t^) ' corresponding to the (alternative) choice of all the master integrals 
with numerator equal to one and with one or two as exponents of the denominators, with at least 
two exponents one. 

9.2 The test case 

As first test of our approach, we have considered in particular detail the case where all masses arc equal, 
mi = ... = fnNd — 1; ^iid all the external lines are on-mass-shell. Denoting the incoming external 
momenta -pi, pi-p2, ■ ■ ■ , PNj,-i -PN^, Pn,,, we choose pf = -1 and {p^ -pj)"^ = -1 for every i and j. 
In the case of box diagrams this corresponds to setting the Mandelstam variables s = t = 1 and u = 2. 
These values of masses and momenta have been chosen because they introduce symmetries which allow 
several consistency checks on equations and results. Similarly to the arbitrary case, in the right half of 
Table 3, for each diagram we list the number n'j, of master integrals with all the denominators, the order 
R' of the difference equation satisfied by the function Ugi{x), the order S' of the differential equation 
satisfied by the corresponding function Vgi{t) and, if present, the number C of non-zero constants needed 
to determine the solution. 

For these particular values of masses and momenta the equations turn out to be, as expected, simpler 
than in the arbitrary case. The homogeneous parts of the equations which presented some similarities 
in the arbitrary case (see the observation 4 of section 9.1) here turn out to be identical; of course the 
nonhomogeneous parts of these equations are quite different. 

The number of constants r]j to find turns out to be greater than or equal to that of the arbitrary case 
(except for diagram 2e); in this connection the test case is more complicated than the arbitrary case. 
With the exception of the diagram 4d6, the calculation of the scalar master integrals requires no more 
than two constants, easily determined using the identities of section 5.5 (which turn out always to provide 
one useful relation involving the constants) and the large-x leading behaviours. One can prove that the 
values of external momenta are always below or at the deformation threshold, with the exception of some 
subdiagrams of diagrams 4f and 4i; in the cases at deformation threshold the values of the constants 
Li, L2 and L3 given in Eq.(94) are needed (for example, for the diagrams 2a, 3a and 4a). 

The characteristic equations of all difference equations have always the solution /U = 1. We list the 
values of the other roots for some diagrams: (Fig. lb: /U = —1/3), (Ic: —1/8), (le: —1/2), (2e: —1/3, 
-1/15), (2il: -1/3, -1/2), (2i2: -1/3), (2i3, 2j2, 211 and 3b2 : 1± (2vl-2 and 3f: -1/3, 1/9), 

(3a: -1/2), (4a: -3/5), (4bl-2: -1/2, -1/4 ± i,/l/8), (4b3: 3/4 ± ^27/32), (4fl: -1/2, -1/3, -1, 

The equivalence of recurrence relations very recently described in [20] may probably throw light on this. 
-'^Wc consider here the master integrals containing one particular combination of denominators. The total number of 
combinations (8) itself grows exponentially with A^^, clearly along with the total number of master integrals with every 
possible combination of denominators. 
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(l±\/3)/2, -0.098 ±0.050i, -0.934, 3.632), (4f2: -1/2, -1/3, -1, (1 ± V3)/2, 0.049 ± 0.099i, -2.607, 
-0.819, 1.756, 12.07), (4il: 1/9, 3/2, -1/2, -1/3, -1) and (4i2: 1/9, 3/2, -1/3, -3). Most of these 
values, written in the form l/(g^ + 1), are due to presence of singularities at k"^ = in the function /(fc^) 
(sec section 5.4). All characteristic equations turn out to have one negative rational solution —1 < /i < 0, 
so that, according to section 6.1, all recurrence relations are unstable; the diagram 2e shows the greatest 
instability {A = 15). One can ask whether the difference equations of all the analyzed diagrams admit 
as solutions convergent factorial series expansions. The answer is negative. In fact the root /i = 1/9 of 
the diagrams 2vl-2 and 3f, the root (1 + •\/3)/2 and the complex roots of diagrams 4fl-2, the roots 1/9 
and 3/2 of diagrams 4il-2 (note, all the diagrams with lines crossed) and the root 3/4 + -\/27/32 of the 
diagram 4b3 satisfy the condition (102) (with f/°'^ = 1) so that the factorial series expansions of the 
solutions never converge. Therefore we are forced to use Laplace's transformation for these diagrams and 
these choices of the line Di, and for all the diagrams which become these diagrams by deleting lines. In 
all the other cases factorial series expansions can be used. 

9.3 Numerical results 

For each diagram shown in Figs. 1, 2, 3 and 4, except for the diagrams 4f and 4i (see section 9.3.4), we 
have calculated the values of the master integrals for £) = 4 — 2e, using the values of masses and momenta 
shown in the previous section. 

Calculations were carried out by using the program SYS described in section 11. The master integrals 
of the diagrams lb, le, 2a, 2d, 2i, 2j, 21, 3a, 3e, 4a and 4c were first calculated using expansions in factorial 
series; these integrals were also recalculated using Laplace's transformation in order to provide important 
checks of the calculations. The master integrals of all the remaining diagrams were calculated using 
Laplace's transformation. Excluding the simplest diagrams, calculations with Laplace's transformation 
turned out faster than calculations with factorial series; this is due to the instabilities of the recurrence 
relations, which become deeper increasing the number of loops, and which force calculations with factorial 
series to be performed with a larger number of digits. To give an idea of the size of calculations, the 
systems of difference equations between the master integrals of the diagrams le, 2d, 2t, 2u, 2v, 3f, 4g, 
and 4h are formed with 44, 28, 245, 304, 362, 81, 139 and 158 equations, respectively; note that in each 
system, in the right-hand side of the equation corresponding to the more complicated mastcir integral 
almost all the other master functions appear. We made no use of the symmetries due to the particular 
values of masses and momenta in order to simplify or reduce the number of the equations, as the aim 
of program SYS is to deal with calculations of multi-scale integrals, lacking in such symmetries; we used 
them only to check the final results. In order to guarantee results with at least 20 digits of precision, 
calculations with factorial series were performed with precision up to 77 digits (depending on the degree 
of instability), while calculations with Laplace's transformation were all performed with 38 digits of 
precision. 

For example, the calculation from scratch of the integral /(2t), Eq.(220), with Laplace's transforma- 
tion, requested about 128 hours of CPU time on a 133 MHz Pentium PC; 16 hours wore used for the 
determination of the systems of difference and differential equations, obtained by solving systems up to 
43000 identities. The solution of the systems (245 equations) yields, as a byproduct, also the values of all 
simpler master integrals, including Eqs.(205)-(208), Eqs.(210)-(214), Eq.(216) and Eq.(219). We stress 
that at this preliminary stage of development we directed our efforts to devise tests and cross checks 
rather than to speed up the program. 

For brevity, for each diagram we list here only the values of the master scalar integrals 



without scalar products and containing all the Nd denominators. As usual, the results have been nor- 
malized with the division by = r(l -|- e) raised to the number of loops of the diagram. Coefficients 
are shown with only 13 digits to save space. Values for e = of all finite integrals have been checked 
by comparing them with numerical values obtained by performing Monte-Carlo integrations over Feyn- 
man parameters, or by performing low dimensional Gaussian integrations on integrands obtained using 
dispersion relations and hyperspherical variables. As additional consistency check we have repeated the 
calculation of some diagrams with different choices of the line Di , and we have checked that the results 
obtained are the same. 




(195) 



43 



9.3.1 Vacuum diagrams 



I{l&)rj^ = -e-i - 1 - e - - - + O(e^) , (196) 

I{lh)T-'^ = -1.5e-2 - 4.5e-^ - 6.984139141966- 18.00878162355e 

- 27.99422356368e2 - 72.003786597996^ - 111.9974983355e^ + 0{e^) , (197) 



I{lc)r-^ = 2e-^ + 7.666666666667e-2 + 17.56"^ + 22.91666666667 

- 661.1105861534e^ - 3685.054779382e* + 0{e^) , (198) 



+ 21.25179105129e- 184.2300051053e^ 



/(ld)r73 = -e-^ - 5.666666666667e"2 _ I5.30161161726e^i 

- 394.1378145809e^ 

- 1375.66943521 le^ - 3466.9749998996e^ + 0{e^) , (199) 



- 46.07511172933 - 148.3508545129e - 394.1378145809e^ 



I{le)T;^ = 2.404113806319e-i - 10.03527847977+ 35.94478903214e 

- 119.1503507802e^ + 379.7433345095e^ - 1183.320931551e'' + O(e^) . (200) 

Eq.(197), first six terms of Eq.(198) and first two terms of Eq.(200) agree with the analytical expressions 
given in [21], [4] and [22], respectively; remaining terms and other results are new. 

9.3.2 Self-energy diagrams 

/(2a)r7i = e"^ + 0.186200635766+ 0.021156303568e+ 0.0017267453536^ 

+ 0.000109897792e^ + 0.000005730593e^ + O(e^) , (201) 

I{2h)T-'^ = -1.5e-2 - 4.25e-^ - 7.375 - 17.22197253479e 

- 29.559207053726^ - 68.87789517038e^ - 118.2464846454e'' + O(e^) , (202) 

I{2c)T-^ = G.5e~2 + 0.68620063576586-^ - 0.6868398873414 
+ 1.4863983919136- 2.9387965877456^ + 5.8710863659586^ 

- 11.736165714496-* + 0{e^) , (203) 

I{2d)T-'^ = 0.9236318265199- 1.2849216718486+ 2.6895076264906^ 

- 5.3383992275116^ + 10.671367369126* + 0(6^) , (204) 



I{2e)T-^ = 2e-^ + 7.333333333333e"2 ^ 16.027777777786-* 

184.141366543 

- 838.23643241786^ - 3647.1021970316^ + 0{e^) , (205) 



+ 21.92956264368+3.6051274751616- 184.14136654316^ 



I{2f)r-^ = -6"^ - 2.6126342869826-2 - 3.9064204906906- * 
+ 0.5840769314959+ 1.764600414536+ 107.00720314356^ 

+ 163.28552937836^ + 1372.2414661896* + 0(6^) , (206) 
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I{2g)r-^ = -(-^ - 5.5g-2 - 15.48413914197e-^ - 45.68793012675 

- 149.1607636537e- 392.298867227e2 

- 1380.1258331676^ - 3455.548194007e'* + O(e^) , (207) 

I{2h)T-^ = - 5.333333333333e-2 - 16e"i - 43.91483126325 

- 154.918028663e- 374.09418533346^ 

- 1436.6727125356^ - 3281.940436319e^ + 0{e^) , (208) 

I{2i)T-^ = 2.4041138063196"^ - 9.763424447585 + 34.99888165588e 

- 116.04204775646^ + 370.04072740696^ - 1153.646312515e^ + O(e^) , (209) 

I{2})r-^ = 0.33333333333336-3 + 0.51953396909916"^ 

+ 0.57536094942696"^ - 2.981838135558+ 12.565962041086 

- 44.853023515386^ + 149.17428117216^ - 477.14408861296^^ + 0{e^) , (210) 

I{2k)T-^ = 0.16666666666676^3 + 0.59310031788296"^ 

+ 0.062348945424026"^ - 1.364667486582+7.0624824278946 

- 26.874199155736^ + 91.916412844176^ - 298.1969436136* + 0{e^) , (211) 

I{21)TJ^ = 0.16666666666676"^ + 0.59310031788296"^ 

+ 0.062348945424026"^ - 1.158877567105+6.2686605834276 

- 24.11937497596^ + 82.98720593436^ - 270.16887601036* + 0(6^) , (212) 

I{2m)r;^ = 0.33333333333336^3 _^ 0.85286730243246^^ 

- 1.7281694115846^1 + 6.070141409747- 19.486513655166 

+ 61.388287566276^ - 190.330269530663 + 583.80453815296^' + 0{e^) , (213) 

/(2n)r73 = 0.33333333333336"3 + 0.85286730243246"^ 

- 1.7281694115846"^ + 6.120359708375- 19.670630424676 

+ 62.001782532356^ - 192.25862531846^ + 589.72125447166* + 0(6^) , (214) 

/(2o)r73 = 0.92363182651996"^ - 2.423491634417+8.3813497100696 

- 26.993621216776^ + 85.100963229996^ - 263.9033186296* + 0(6^) , (215) 

/(2p)r73 = 0.92363182651996-1 - 2.116169718457+6.9295446852596 

- 21.503278377386^ + 66.3221338040163 - 202.8870257176* + 0{e^) , (216) 

/(2q)r73 = 1.326448208272- 5.1966481369656+ 18.377583878046^ 

- 60.411915036616^ + 191.59639416* + 0{e^) , (217) 

/(2r)r73 = 1.341399241447 - 5.1977529558966+ 18.387046564076^ 

- 60.423352130163 + 191.6140096256* + 0(6^) , (218) 
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/(2s)r7^ = 2.002500041105- 8.162562835907e+ 29.46716463085e2 

- 98.13080591819e3 + 313.871881187e'' + O(e^) , (219) 

I{2t)rj^ = 0.2796089232826- 0.1380294113932e+ 0.3194688268113e^ 

- 0.4399664109267e3 + 0.6650515012166e'^ + O(e^) , (220) 

J(2u)r7^ = 0.1826272375392- 0.06746690965803e+ 0.18654624206236^ 

- 0.2498713405447e^ + 0.3796187113121e'' + O(e^) , (221) 

7(2v)r7^ = 0.1480133039584- 0.009263002043238e+ 0.10533085373976^ 

- 0.1224292041846e^ + 0.1898480457555e^ + O(e^) . (222) 

Eq.(202), first four terms of Eq.(203), first term of Eq.(204) and first six terms of Eq.(208) agree with the 
analytical expressions given in [23], [24], [25] and [4, 5], respectively; remaining terms and other results 
are new. 

9.3.3 Vertex diagrams 

/(3a)r7^ = 0.671253105748+ 0.1998957762816e+ G.03189366853371e2 

+ 0.003532937320333e^ + 0.0003018185047825e'' + O(e^) , (223) 

I{3h)T-'^ = 0.5e-^ + 0.68620063576586"^ - 0.5916667014024 
+ 1.3561965331146- 2.6691121188146^ + 5.3366513585166^ 

- 10.668662837416^ + 0{e^) , (224) 

I{Sc)r-'^ = 0.6712531057486"^ - 0.08774519609257+0.72623756269476 

- 1.321129485876^ + 2.6674314693766^ - 5.3326753370916^' + 0(6^) , (225) 

I{3d)Tj^ = 0.937139527315- 1.271849687086+ 2.6918504750662 

- 5.3369321349616^ + 10.671003429346^' + 0{e^) , (226) 

/(3e)r7=^ = 0.2711563494022+ 0.18339410775146+ 0.053751010587696^ 

+ 0.014461033684196^ + 0.0007461873722766^ + 0{e^) , (227) 

7(3f)r72 = 0.173896742268+ 0.18166648769626+ 0.044408991818326^ 

+ 0.022315473857856^ - 0.0030798104797976* + 0(6^) . (228) 

9.3.4 Box diagrams 

7(4a)r7^ = 0.3455029252972+ 0.47310083188186+ 0.15194595375436^ 

+ 0.02751792845546^ + 0.003484921775196* + 0{e^) , (229) 

7(4b)r72 = 0.6712531057486"^ - 0.06425178040942 + 0.7393966927045e 

- 1.3171876991126^ + 2.6681073433996^ - 5.3325008824596* + 0{e^) , (230) 
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/(4c)r7^ = 0.9509235623171- 1.258189955235e+2.694588643167e2 

- 5.335347124508e3 + 10.67067120761e'' + O(e^) , (231) 



I{M)T-'^ = 0.276209225359+ 0.1937422320842e+ 0.06034849310181e^ 

+ 0.016408285886816^ + 0.0013016425167656^^ + O(e^) , (232) 

J(4e)r7=^ = 0.3455029252972e-^ + 0.4347670080988+ 0.17885718095363e 

+ 0.05005382113385e^ + 0.006936292250698e^ + 0.002511559375421e^ + O(e^) , (233) 

I{Ag)T-'^ = 0.1723367907503+ 0.2679578491711e+0.13552112755141e2 

+ 0.04468531532833e^ + 0.008430602827459e^ + O(e^) , (234) 

7(4h)r72 = 0.1036407209893+ 0.2142416932987e+0.14046068671363e2 

+ 0.04437197236388e3 + O(e^) . (235) 

In the case of the diagrams 4f and 4i (from which the diagram 4f is derived with the deletion of a Une) 
the equations (27), obtained after the transformation of the subsystems of equations into triangular form, 
have orders relatively high (> 10) and expressions of size greater than the limit of the computer used, 
so that our program failed to work out them (the coefficients of the equations would be polynomials in 
two variables of degree ~ 100); we were able to find only the characteristic and indicial equations. The 
calculation of these integrals without the transformation into triangular form will be considered in the 
next paper. 

9.4 Some examples of equations 

Here we discuss the equations satisfied by the functions Ig{x) = Ugi{x) for some simple diagrams with all 
lines topologically equivalent, using the values of masses and momenta of section 9.2. For the diagrams 2a 
and lb the difference equations are 

Wi72a(a;) = -liD- 2)/ia(x) , (236) 
WiJib(a;) = -{D - 2)/ia(l)/ia(a:) , (237) 

where Wi is the operator 

Wi7(a;) = -3xl{x + l) + {2x-D + l)I{x) + {x - D + l)I{x - 1) , (238) 

and Iiii{x) is the integral (103) with m = 1. The identity of the homogeneous parts of the equations is an 
example of the similarities between equations described in section 9.1. The solutions of Eqs.(236)-(237) 
can be written as 

l2.{x) = mJi'J' (x) + Iii"{x) , (239) 
/i,(x) = mi^liL'^ix) + 2/ia(l)/f/(x) . (240) 

Numerical calculation of J|^'^(l) and /2a^(-'-) '^^^ described in section 7.2.4. The constant rjih may be 
found by comparing the large-x behaviours of the solutions /2a*^(^) ~ a;~^/^+^/^, lH^ix) cx x^^^^ (see 
Eqs.(126)-(127)) with 



(241) 



(see Eq.(71)). One finds r]ih = 0. Therefore /ib(l) = 2/ia(l)-f2a^(l)i that the two-loop integral 

/ib(l) factorizes into a product of a vacuum one-loop integral and a part of a one-loop self-energy integral. 
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Now let us consider the diagrams of Figs. 2b and Ic. Both diagrams have only one master integral 
with all the denominators, the scalar integral (compared with the four master integrals of the arbitrary 
case); the corresponding master functions l2h{x) and /ic(a;) satisfy the difference equations 

W2/2b(a:) = {D- 2)2/i,(l)/i4x) , (242) 
W2/ic(x) = §(£> - 2fll{l)h,{x) , (243) 

where W2 is the operator 

W2/(a;) = -8x{x -D + 2)I{x + 1) + {ix^ + (13 - 1015)3; 

+ (3£> - A){D - l))l{x) + {x-D + l){x- 'iD/2 + 2)/(x - 1) . (244) 

Note again the identity of the homogeneous parts of the equations. The characteristic equation has the 
solutions 1 and —1/8. The index associated to /x = 1 is K = —D/2; therefore the condition (73) is 
satisfied and the corresponding solution of the homogeneous equation I^{x) contributes to /2b (a;) and 
I\c{x). The solutions of Eqs.(242)-(243) can be written as 

h^{x) = V2^lTix) + I^^"{x) , (245) 
hc{x) = mJ^^Fix) + lh.{l)I^^"{x) . (246) 

The constants may be obtained by comparing the large- a; behaviours I^^(x)~ and /^^(x) oc 

with J2b(x) ^ x-^/2/2a(l) and hcix) « h^{l). One finds 7?2b = hail) and Vic = /ib(l)- 

We consider the calculation of 1^*^(1) and /^^(l) with factorial series. The recurrence relation (244) is 
unstable with A = 8. We fix a precision of the coefficients of the powers of e of the results of -E = 13 digits, 
with = 7 terms of the expansions in e; following section 6 we choose Xmax = 25, and we perform the 
calculations with C = 38 digits. Expansions in factorial series of I^^{x) and I^^{x) converge for x = 25 
in about 1800 terms. Values for x = 1 of the solutions are calculated using repeatedly the recurrence 
relations (242). Two terms of the expansions in e are lost going beyond the abscissa of convergence A = 3, 
so that we must retain the first rig = 9 terms of the expansions. One obtains 

V2hl2h° {1)1^7'^ = 0.09188814923697e~3 + 0.194632539439e-2 

- 0.045490472375e-i + 6.50912255436 + 10.43240978278e 

+ 76.7023111407e2 + 118.6149739413e^ + 732.0021 187015e^ + 0{e^) , (247) 

7^^(l)r72 = -0.09188814923697e-3 - 1.694632539439e-2 

- 4.204509527625e-^ - 13.88412255436 - 27.65438231756e 

- 106.2615181944e^ - 187.4928691117e^ - 850.2486033469e'' + O(e^). (248) 

Summing the results one finds Eq.(202); note the cancellation of the term, and the partial cancellations 

of digits in the e^^, e^^ and constant terms. Considering now the solution of Eq.(242) with Laplace's 
transformation, there are three master transformed functions: V2h(t), wi(t) and W2{t) defined by 







Clearly wi{t) ^ W2{t) and Ki{x) ^ K2{x) as the insertion of x breaks the symmetry kx ^ k2- The 
function V2h{t) satisfies the differential equation 

i^(l -t){l + 8i)w2b + *(8(1 - -D)t^ + W{D - 2)t + (5£)/2 - 6))u2b 

+ (3D - 8){D - 3)(i + l/2)t;2b = [D - 2ftl2a{l)v2a , (251) 
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where f2a is v,j of Eq.(184) with m = 1. i^i(l) and i^2(l) arc not master integrals, as Ki{\) = K2{\) = 
— ^/2b(l); therefore Ki[x) and K2{x) are not master functions and satisfy difference equations of order 
zero, for example 

(x + 2 - D){~iKi{x + 1) - /2b(x + 1) + /2b(a;)) = (2 - D)h^{x + l)/2a(l) • (252) 
On the contrary, Wi{t) and W2{t) satisfy first order differential equations, for example 

- &t{tw[ + {D- 2)wi) + 2t{l - i)w^b + 2((2 -D)t + D- 3)v2h = 2(2? - 2)72a(l)to2a , (253) 
and therefore they are master functions. 



10 Differential equations in masses and momenta 

As shown in the one-loop example, the calculation of the constants rjj , the factors multiplying the solutions 
of the liomogcncoTis equation, presents different degrees of difficulty according to the values of masses 
and momenta. For values below the deformation threshold, these factors (when non-zero) are easily 
expressed in form of integrals with one loop missing, but, above the deformation threshold or when some 
masses arc zero, as discussed in section 5.4, the calculation of the factors becomes more complicated. 
An alternative to this calculation is to calculate the master integrals for values of masses and momenta 
such that the calculation of the factors is simple, and subsequently to build and integrate a system of 
differential equations in the masses and momenta in order to reconstruct the integrals for the requested 
values of masses and momenta. The approaches to calculation of diagrams based on construction of 
differential equations in masses or momenta were introduced in [26, 27]; up to now they have been used 
in the analysis of particular diagrams (e.g. massless box diagrams in [9]) with a small number of master 
integrals by building and solving small systems of differential equations. Here we describe how to build 
and solve systems of differential equations in masses and momenta in way very similar to that of the 
differential equations obtained by Laplace's transformation. In this way even large systems containing 
hundreds of equations can be built and solved, obtaining results expanded at will in e. 



10.1 Differentiation of master integrals 

A generic master integral is a function B(m^,pi ■ pj) of scalar products and masses. We parametrize 
the trajectory from the initial point to the final point in the space of scalar products and masses with a 
parameter y; the scalar products and masses of the points of the trajectory will be in general functions 
of y. In analogy with the di£Ferc;ntial c;quations obtained by using the Laplace's transformation, we 
introduce the parametrization in such a way that the initial point corresponds to y = 1 and the final 
point corresponds to y = 0. Derivatives with respect to y of master integrals will be written as derivatives 
with respect to scalar products and square masses: 



dB dpi dB ^ y d{p, ■ p,) dB ^ ^ dmi dB 

dy 2-, dy dpi f-' dy d{pi-pj) fri dy dmi' 



Differentiation with respect to masses of the integral is trivial; differentiation with respect to scalar 
products is less immediate. Let us consider a diagram with Np independent external momenta p,; the 
scalar products pa ■ Pp are Npr = Np{Np + l)/2. The derivative with respect to pa ■ pp has the general 
form 

d d 
d{Pa-P0) dpj 

for (a, /?) equal to each one of the Npr possible different pairs of integer. The coefficients aya/J can be 
determined by imposing the Npr^ conditions 

djp-y ■ Ps) ^ h^tP'r-P5=Pa-P0 , .25g^ 
9{Pa • P/3) jo otherwise . 
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But the number of the coefficients aijap is Np^Np,. which is greater than Npr'^; therefore some of these 
coefficients are not constrained by these conditions and may be chosen arbitrarily. We conveniently 
rewrite Eq.(255) as 

where the coefficients aijap to be foimd are Np^'^ and there are Np{Np — l)/2 arbitrary constants Adij. 
In order to calculate the coefficients it is convenient, fixed a and /?, to apply Eq.(257) to a generic scalar 
product Pk-Pi, and to construct the tensor Aiji-i formed by the factors multiplying ttya/J in the obtained 
expression. Now, considering the various pairs of indices {i. j} as one single index ij, A becomes a square 
matrix Aij^ki of dimension Npr- Inverting it we find the desired coefficients 

As an example, for a vertex diagram Np = 2, and the matrix A is 

/ 2pl 2Mi2Pi-p2 \ 
A = \pi-p2 pI + Mi2pI Pi • P2 , 
V 2pi-p2 2pl J 

where the indices 1,2,3 correspond to ij = 11, 12, 22. 

10.2 Construction of the system of differential equations 

The algorithm used for the generation and solution of the system of identities is the following: 

Algorithm 4 Follow the algorithm 1 with the following modifications: 

1. New identities which transform the derivatives with respect to y of master integrals into combinations 
of integrals obtained using Eq.(254) are added to the system. 

2. Master integrals and their derivatives have a priority of extraction lower than other integrals. 

3. Add the new entry "the greatest derivative" to the list of priorities after the entry 9(b)vi. 

Analogously to sections 2.5 and 3.2, the previous algorithm allows one to determine the list of the 
master integrals Bn{y), and to obtain a system of differential equations between them. Applying a 
procedure of transformation into triangular form of the subsystems of equations between master integrals 
with the same denominators, the whole system of differential equations takes the triangular form. The 
arbitrary constants Mij turn out to be particularly useful as a check, because the expression of derivatives 
(254) contains Mij (through Eq.(257)), while the differential equations obtained from the system of 
identities are obviously independent of them. 

The integration of the system of differential equations requires the knowledge of the values of the 
master integrals and of some derivatives in the initial point y = 1. The values of the master integrals 
are obtained by solving a system of difference equations. The derivatives of the master integrals in the 
initial point are expressed in terms of the master integrals of the system of difference equations using the 
identities provided by the algorithm 1. 

Note that the master integrals of the system of difference equations may be different in number and 
structure from the master integrals of the system of differential equations. The system of differential 
equation is solved with the same procedure used for the Laplace's transformation method. We stress that 
mobile apparent singular points, as described in section 8.5.2, may appear in the system of differential 
equations. 

10.3 Expansion in e and infrared divergences 

All quantities, coefficients of differential equations, solutions, T hmctions, etc. are expanded in e = 
(4 — D)/2. Solution of differential equations using power series in y with coefficients expanded in e does 
not cause difficulties, unless we deal with IR divergences. If a master integral B{y, e) is IR finite in the 
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Figure 5: Diagrams with massless lines (dashed). 

initial point y = 1, but is IR divergent in the final point j/ = 0, it is not possible to obtain the correct 
dimensionally regularized value, which contains an additional IR pole in e, by integrating the differential 
equation from the initial point, where the integral is IR finite: the coefficients of the expansion in e 
of B{y,e) tend to infinity as y ^ 0, so that limy^Q B{y,€) ^ B{0,e). In this unfortunate case the 
only remedy is to substitute the master integral IR divergent with another integral (or combination of 
integrals) IR finite. Our preference for master integrals with scalar products in the numerator and all 
exponents 1 in the denominator rather than integrals with exponents greater than 1 in the denominator 
was devised to limit the appearance of master integrals IR divergent. 

In the case where B{0, e) is IR finite, but some derivatives with respect to y of B{y, e) arc IR divergent 
in 2/ = 0, the problem is overcome by integrating the differential equation up to a small value y ~ 0, 
suitable chosen so that the value of the master integral will be calculated with the requested precision. 



10.4 Examples of calculation of integrals with zero masses 

For illustrative purpose, we consider here only two simple diagrams; applications to more complicated 
diagrams will be shown in future papers. Let us consider the diagram of Fig. 5a, and the integral 

H^) = I 772 .^ ,J1^^ , ' = -1 , (260) 

^ ' J {kf + 1)^ {p-ki- k2y ' ' ^ ' 

where [dk] = [d^ki] [d^k2\- It satisfies the homogeneous difference equation 

x{x -2D + 5)I{x + 1) - (a; - 3D/2 + 3){x -D + 2)I{x) = ; (261) 

the solution of this equation is 

The calculation of the constant rj is not simple; in fact, the value of — = 1 is above the deformation 
threshold = 0, so that the large- a; behaviour of I{x) contains an additional contribution due to the 
turning point of the integration path; moreover, the identity (100) is useless because 7(0) = 0. 
We replace the zero mass with a square mass y in the denominators 

Di = kf + l,D2 = kl + y,D3 = {p-ki- kif + y , (263) 

then we build the system of identities following the algorithm 4. The master integrals are Ii{y) = 
J[dk]/DiD3, h{y) = J[dk]/DiD2, hiy) = /[dfel/DaDg, hiy) = J[dk]/DiD2Ds and h{y) = J[dk]{p- 
k2)/ D1D2D3, and satisfy the system of differential equations 

yl[ = {D- 2)/i , 2yl'^ = [D - 2)h , 2y4 = (D - 2)h , 
8y(y - = 4 {{AD - I2,)y + 7 - 2D) ^ - 2(3i? - S){D - 3)h + {D - 2f{h - h - h) , 
6{D - 2)h = -My - 1)^4 + 2{{D-3)y + 5- 2D) h + {D- 2)(2/2 - /i - h) . (^64) 
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The values of the integrals in the initial point y = I can be expressed using the values of integrals with 
masses equal to one given in section 9.3: /i(l) = hil) = hil) = /^(2a), h{l) = /(2b), 15(1) = -/(2b)/3. 
The singular points of the system are y = 1 and y = 0. Considering the equation for I4, the exponents 
in these points are 

p^yr'^ = 0, 2 - 2e , = 0, 3/2 - 2e . (265) 

In the initial point the integral /4(1) and its derivatives are IR finite, so that the solution must be regular; 
therefore the singular solution associated with p^~^^ does not contribute, and the second-order equation 

for 74 needs only 14(1) as initial condition. In the final point the integral /4(0) is IR finite but ^4(0) is 

IR divergent, so that the singular solution associated with p'^~^^ must contribute to l4{y). 

The integration of the system is divided into two parts: from y = ltoy=l/2 expanding in y = 1, 
and from y = 1/2 to the value y = <C 1 expanding in y = 0. The presence of the cutoff A is needed 
because the solution is not regular in the origin. Perusing the numerical results, the effect of the cutoff 
on the coefficient of of the expansion in e of the final value turns out to be approximately proportional 
to A^ log™^* A, where m is some small integer. A value A = lO"""^^ suffices to obtain the first coefiicients 
with 20 exact digits. The normalized result is 

h{0)r-'^ = -0.5e-2 - 1.25e-i - 4.6648681336964528729 
- 9.595397946879509324e- 26.045799878017610383e2 

- 49.5019341875628515466^ - 120.38235865133474218e'^ + O(e^) . (266) 

Another more complicated example is the integral of Fig. 5b 

[d^h] [d^fes] Kh] 



{kf + 1) ((fci - fc2)2 + 1) ((fe - ks^ + 1) kj ki (fci - ksY 



(267) 



As before, we give a square mass y to masslcss denominators. In this case there are 43 master integrals 
from 2 to 6 denominators, of which only 1 with 6 denominators. The system of differential equations is 
too long to be shown here. The singular points are y = —1, 0, 1/9, 1/4, (3± \/5)/2, 1,4, 9; the points 1/9, 
1/4, and (3 — \/5)/2 on the interval [0, 1] turn out to be regular points of the solutions. The effect of the 
cutoff is different from the previous case, being proportional to Alog™^* A. A value A = lO"'^" suffices to 
obtain the first coefiicients with 20 exact digits. The result is 

LT-^ = 2.4041138063191885708e-i - 3.0270094939876520198 
+ 22.804522068631748454e- 53.102275435449702689e^ 

+ 201.843339942193966946^ - 577.74024368094326834e^ + O(e^) . (268) 

The first two terms agree with the results [28]; subsequent terms are new. 

These calculations were carried out by using the program SYS; the calculation from scratch of the 
integral L and all the other master integrals, including the calculation of the master integrals in the initial 
point, required about 1.5 hours on a 133 MHz Pentium PC. 



11 The program SYS 

Here we report in brief some information concerning the program SYS used to calculate the values of the 
integrals of the sections 9.2 and 10.4. Further details will be given in a separate paper. 

• C program, about 23000 lines, 1Mb executable. 

• The program SYS allows one to calculate the value of integrals in almost completely automatic way; 
only needed input is a description of the diagram, the constants Ui and bi of section 2.6 which limit 
the identities generated, the number of dimensions Dq about which to expand the integrals, and 
the number of terms of the expansions in D — Dq. 

• The program contains a simplified algebraic manipulator, used for the solution of the systems of 
identities. 
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• CocfRcicnts of the integrals in the identities can be unlimited precision integers, rationals, ratios of 
polynomials in one and two variables (for example D and x) with integer coefficients. At present 
the values of square masses and products of external momenta must be rational numbers. 

• Efficient management of systems of identities of size up to the limit of disk space (tested up to half 
million of identities). 

• Numerical solution of systems of difference/differential equations up to 500 equations. 

• Numerical variables used in the solution are arbitrary precision floating point complex numbers and 

truncated series in e with this kind of coefficients. 

• Arithmetic libraries which deal with operations on integers, polynomials, rationals, floating point 
numbers and truncated series in e were written on purpose. 

Two versions of the program exist, one using factorial series and one using Laplace's transformation. The 
Laplace's transformation version is much more complicated than the factorial series version but it turns 
to be faster in many cases: both systems of difference and Laplace-transformed differential equations are 
generated, then the system of differential equations is solved, and the functions U (x) are obtained by 
integrating over t the solutions v{t), also checking that U{x) are solutions of the system of difference 
equations. Both versions of program were used for calculating the single-scale integrals of section 9.2. 

12 Conclusions 

Most part of this paper has been devoted to the description of a new method of calculation of master 
integrals, based on numerical solution of systems of difference equations, obtained by solving systems 
of identities obtained by integration-by-parts. Important features of the method are the applicability 
to arbitrary diagrams, inherited from integration- by-parts method, and the ability to obtain very high 
precision results (even 100-200 digits) expanded at will in e, due to fact that the calculation of integrals 
is reduced to sums of convergent factorial or power series in one variable. 

We also have described two important complements to this method: an algorithm for the reduction of 
generic Feynman integrals to master integrals, and a procedure for construction and numerical solution 
of differential equations is masses and momenta; among other things, at present the latter is needed 
to calculate generic master integrals with massless denominators or with external momenta 'deep' in 
the non-euclidean region (over the deformation thresholds), at least until a working automated general 
procedure for the calculation of the arbitrary constants of difference equations in these cases will be 
found. 

The implementation of these methods and algorithms in the program SYS turned out to be essential 
to test and prove the validity of the approach. 

In this first exploratory work, mostly devoted to test our approach, we have focused our attention 
to the calculation of single-scale master integrals, in particular vacuum and self-energy diagrams up to 
three loops and vertex and box diagrams up to two loops. The calculated values of these master integrals 
are useful, as they may appear expanding Feynman integrals with respect to ratio of different scale 
parameters. 

But the final targets of our approach are the calculation of four-loop g-2 contribution in QED, and the 
calculation of multi-scale multi-loop master integrals especially in cases where there arc no hierarchies 
between scales, where asymptotic expansions in large masses or momenta seem to be not useful. No 
insurmountable difficulty seems to exist for applying our approach to these problems. Clearly, that 
means to deal with a larger number of master integrals, or with more complicated equations, and that 
will imply modifications or improvements of various parts of the algorithms implemented in the program 
SYS which, at this stage of development, is far to be optimal. The experience gained by performing the 
calculations of this work has given many suggestions on the changes which should be made and which 
will be discussed in future papers. 
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Figure Captions 



Figure 1: Vacuum diagrams up to three loops. 
Figure 2: Self-energy diagrams up to three loops. 
Figure 3: Vertex diagrams up to two loops. 
Figure 4: Box diagrams up to two loops. 
Figure 5: Diagrams with massless lines (dashed). 
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Table 1: Values of J{x), J5'^(x) and I^"{x) 
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Table 2: Dependence of the finite part of 7(1) on x. 
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Table 3: Number of master integrals and orders of the equations for the diagrams. 
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